{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import math\n",
    "from numpy.random import randn\n",
    "from scipy.linalg import inv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kalman Filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kalman filter is a bayesian filter used to approximate system state using given measurements and process model. Below I will implement Kalman filter for multivariable system and demonstrate its results. Next step is to improve Kalman filter by implementing dynamic calculation of process noise during filters operation in order to improve its precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Research Process Description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to show results of the research, following plan was developed:\n",
    "1. Develop a model with several variables, such as train movement.\n",
    "2. Create imprecise measurements of the model, in order to use filter on\n",
    "3. Show performance of a classic Kalman Filter on the model and measurements\n",
    "4. Design Kalman filter with uses an agent to calculate process and sensor error as the filtering is in process\n",
    "5. Show differences in performance of the filters - accuracy, speed, memory usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Under Study"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example we will use train for which we have a location sensor. We will estimate speed and location as variables for our model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kalman Filter Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![kalman](img/kalman-2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In principle kalman filter has two stages prediction and update. During prediction we compute _prior_ which is the estimation of the next state of the system. During update step we incorporate measurements and select value of the system somewhere between prediction and measurement based on _Kalman gain_.\n",
    "\n",
    "Below we will see how these stages are described and how we implement them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kalman Filter Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initialization**\n",
    "\n",
    "    1. Initialize the state of the filter\n",
    "    2. Initialize initial state of the system\n",
    "    \n",
    "**Predict**\n",
    "\n",
    "    1. Use process model to predict state at the next time step\n",
    "    2. Adjust belief to account for the uncertainty in prediction    \n",
    "**Update**\n",
    "\n",
    "    1. Get a measurement and associated belief about its accuracy\n",
    "    2. Compute residual between estimated state and measurement\n",
    "    3. Compute scaling factor based on whether the measurement\n",
    "    or prediction is more accurate\n",
    "    4. set state between the prediction and measurement based \n",
    "    on scaling factor\n",
    "    5. update belief in the state based on how certain we are \n",
    "    in the measurement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variables in Kalman Filter\n",
    "\n",
    "Now let's describe all variables we need for a Kalman filter:\n",
    "\n",
    "$\\mathbf{x}$ - state vector, when we talk about state, we always assume mean value, since we are never sure.\n",
    "\n",
    "$\\mathbf{P}$ - covariance matrix of our state.\n",
    "\n",
    "$\\mathbf{F}$ - state transition function.\n",
    "\n",
    "$\\mathbf{Q}$ - process covariance, this is uncertainty you have about your model calculations.\n",
    "\n",
    "$\\mathbf{B}$ - input function, transformation of an input vector corresponding to a model.\n",
    "\n",
    "$\\mathbf{u}$ - input vector.\n",
    "\n",
    "$\\mathbf{H}$ - measurement function.\n",
    "\n",
    "$\\mathbf{z}$ - measurement mean, basically value we get from the sensor.\n",
    "\n",
    "$\\mathbf{R}$ - measurement covariance matrix, how our measurement can deviate from the real value\n",
    "\n",
    "$\\mathbf{K}$ - Kalman gain, or how much we trust either measurement or prediction\n",
    "\n",
    "$\\mathbf{y}$ - residual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equations\n",
    "\n",
    "Kalman filter is described by the following equations:\n",
    "\n",
    "### Prediction\n",
    "\n",
    "$\\bar{\\mathbf x} = \\mathbf{Fx} + \\mathbf{Bu}$\n",
    "\n",
    "$\\bar{\\mathbf P} = \\mathbf{FPF}^\\mathsf T + \\mathbf Q$\n",
    "\n",
    "### Update\n",
    "\n",
    "$\\mathbf y = \\mathbf z - \\mathbf{H\\bar x}$\n",
    "\n",
    "$\\mathbf K = \\mathbf{\\bar{P}H}^\\mathsf T (\\mathbf{H\\bar{P}H}^\\mathsf T + \\mathbf R)^{-1}$\n",
    "\n",
    "$ \\mathbf x = \\bar{\\mathbf x} + \\mathbf{Ky}$\n",
    "\n",
    "$\\mathbf P = (\\mathbf I - \\mathbf{KH})\\mathbf{\\bar{P}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation for testing\n",
    "\n",
    "For testing purposes we need to generate simulation. In our case it is movement of a train. We need to generate measurement with some deviation, as well as our prediction with some process deviation. We will have length of the simulation as a parameter and make samples every second. System will have initial position in form of a vector - first value is initial coordinate, second is initial speed. Each iteration we will change coordinate by speed + process noise. Also when we take measurements, sensor noise will be added. Function will return two vectors - vector of predicted coordinates and vector of measured coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def generate_simulation(initial_state, process_var, sensor_var, length=1):\n",
    "    x = initial_state[0]\n",
    "    v = initial_state[1]\n",
    "    \n",
    "    process_std = math.sqrt(process_var)\n",
    "    sensor_std = math.sqrt(sensor_var)\n",
    "    \n",
    "    coordinates, measurements = [], []\n",
    "    \n",
    "    for _ in range(length):\n",
    "        x += v + randn() * process_std # compute new coordinate\n",
    "        coordinates.append(np.array([[x], [v]]))\n",
    "        measurements.append(np.array([[x + randn() * sensor_std], [v + randn() * sensor_std]])) #compute sensor readings\n",
    "    \n",
    "    return np.array(coordinates), np.array(measurements)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing Simulation\n",
    "\n",
    "Now let's see how results of our simulation look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd8zdcbwPHPyRJ7i71H7ZHYm6oWtdvyQ81qqaKtVlul1VaVqtFSrU21RtUoqiiCVuy99wgRxMiSfX5/nEsSgojce5Pc5/165XXv/d7vOE8u98kZ33OU1hohhBDiQU72LoAQQoiUSRKEEEKIBEmCEEIIkSBJEEIIIRIkCUIIIUSCJEEIIYRIkCQIIYQQCZIEIYQQIkGSIIQQQiTIxd4FeBa5cuXSRYsWTdKxISEhZMyYMXkLlAo4YtyOGDM4ZtyOGDM8fdx79uy5obXO/aT9UnWCKFq0KLt3707Ssd7e3jRq1Ch5C5QKOGLcjhgzOGbcjhgzPH3cSqkLidlPmpiEEEIkSBKEEEKIBEmCEEIIkSBJEEIIIRIkCUIIIUSCJEEIIYRIkCQIIYQQCbJaglBKzVJKXVNKHY6zLYdSar1S6pTlMbtlu1JKfa+UOq2UOqiUqmatcgkhRLK7eBH8/GxyqcC7EYz5aj4Xdx+x+rWsWYOYA7z4wLaPgA1a61LABstrgJeAUpafvsBUK5ZLCCGSh9YweTKULg2dOln1UpHRMczddp5GYzYyNTg7m733W/V6YMUEobXeAtx8YHMbYK7l+VygbZzt87SxHcimlMpnrbIJIcQzu34dWreGd96BjBlh2zYIDX3288bEPLTppH8QzSdu4bM/j1DaKYxVcwbR7cXKz36tJ7B1H4SH1toPwPKYx7K9AHApzn6+lm1CCJHy+PpCjRqwbh1MmgS//AJRUbB9+7Odd/NmyJwZDt9vmWfH2QA6Tt1GUFgUM173YoHfWircvQ7lyj1jEE+WUuZiUgls0wnuqFRfTDMUHh4eeHt7J+mCwcHBST42NXPEuB0xZnDMuG0Rs+utW1QZPJh0N25wYMIEgsqVwzkkhHpOTlyYN4/zTkn/u7vM2LHkCw3l6pAhHB06lF1Xo5l+MJzcGRTvVXXF5doxQtavJ6J0aQ5u2XL/OKvFrbW22g9QFDgc5/UJIJ/leT7ghOX5z0DnhPZ73I+np6dOqk2bNiX52NTMEeN2xJi1dsy4rR7zrVtaV6midfr0Wm/ZEv89T0+tGzdO3Hn+/lvr5s21DguL3RYZqY+Urqrbd/1W13h7ri7+0SpdZOgq3eHH//StkHCzT2Cg1k5OWg8fHu90Txs3sFsn4jvc1k1MfwLdLc+7AyvibH/dMpqpFnBHW5qihBAixXjtNThyBJYuhfr147/XoAH4+EBExJPPM3YsrF0Ly5bd3xTk/S/9GvXjQsGSNDy7l376IqPbV2R+n5pky+Bmdtq92/RR1K6djEE9mjWHuS4AfIAySilfpVRv4BugmVLqFNDM8hrgL+AscBqYDvS3VrmEECJJAgNNn8PQofDigwM0MQkiLMx8iT/OpUuwaZN5PtUM2NRa89HaM/hm9WBqt+qMzXSFIT99TOey2XF3dY491sfHPNasmQwBPZk1RzF11lrn01q7aq0Laq1naq0DtNZNtdalLI83LftqrfXbWusSWuuKWuukLfIghBDWsm+feaxTJ+H369Uzj3H6BhL0669meOxbb5l9jx5l/vYLrHbOy/vXd1G9XAH44AO4cwdmzox/rI8PlCkDOXI8WyyJJHdSCyFEYuzZYx6rPeI+3ly5oHz5+Ali8WL45BOTEMA8zpvH1cbN+fnlfoxq+gbvzfqXL1ceoeHZ3bxVr5jZr3p1UyOZMMGMjrp37PbtNmteAkkQQgiROHv2QIEC4OHx6H0aNIB//zVf6rt3Q9euMHo0/PyzeX/vXo7duEub2v0ZveUS86u1ZMfddNSMusn4vybi1Prl2HN98IG5Q9vSDMWZM3Djhk0TREoZ5iqEECnb3r3g6Xn/ZXB4FCeuBnL5dhh+t+9yKzSS4GLNCGmQnnwz/qHrD5+QP29eKFUK3nsPGjZk26+reLPLGDKmT89fvWpS7swBaNAWnJ1NE1XuOMtEt2wJzZubGkjr1rH9D5IghBAiBQkKghMnoHNnAPzu3KX9j9vwuxN2fxc3FycyuaQnQ8Fy+J2NYFqLz2hZJAMNyuXj8ugJnP/6D1bmrUKx6CDmDKhH/mzpIV890yx15Ai0axf/mkqZmkf58qa/olgxcxOdDW6Qu0cShBDCcfj7w6efwrhxkDVr4o/bv9/0AVSrRnB4FL3m7CYoLIofu1SjZJ5M5MvqTmZ3V7NvyZJcuh7EnPfGsei6OyvWXQDP9uQNukHzk9sY1b0uWbOlN/sqBYMGmek6HkwQAEWKmCaqgQMhfXrTQe7s/PB+ViIJQgjhOKZNgxkzTDNNr16JP87SQR1VtRpv/7qXk/5BzOpRnYalcz+879tvU+j4cYaP6Mrg8Cj8A8MomD0D7h8OgQPLoNWX8ffv08ckh1y5Er52//6wYIFpYrJh8xJIJ7UQwpEsWmQeV69++D2tzZTd69bBjz/CtWux7+3dS3Chonz43zU2n7zOV20rJJwcAN591zQNKUVmd1dK5sls7mUYP950NLu5xd9fqUcnBzA1hhkzoGhRaNXqqcJ9VlKDEEI4hsOHTVt/jhywfr254/nel/W1a6aT+NSp2P19fOCXX4iKjmHhDVcmvjqWG3svM7BpKTrXKPz011cq6c1D5crBuXNJO/YZSA1CCOEYFi4EJycYM8Z0Om/dGvverFkmOXz7rbnL+e23ifltAav+3s2LEzbzaaX2FHeNYsXbdXmvWWn7xWBjUoMQQqR9WpsE0bSpGYk0YIBpZmra1MxtNH06NGwIQ4YQFhnNWicPpoSX56S3PyUzOfHTslE0H/sRqlA2e0diU5IghBBp3549pv3/k0/M4j6NGpkEMX48bNyIPnuWLR+OZsWi/aw76k9weBQlc+Tkhz/H0qJJRZxP+oCX5xMvk9ZIghBCpH2LFoGra+xQ0pYtzdDR06fxmzmfYZ2/YuO5jGRx96dFxby0qpSfullicJ7cE6ZsMXdP589v3xjsQBKEECJti4kxCaJ5c8ie3Wxr2RI9cCCLZq5mVN4WRLq582nLsnSrXYR0LnE6kgcMMFNzV6tmOpkdjHRSCyHSNh8fM8V2p073N+lixfii44d8pEtS/uoZ1nYoRp/6xeMnB4AhQyBbNtM/4YCkBiGESNvmzYMMGcx8RkB0jObT5YdYUKIBvXYt59OokzjV+DjhY3PnhvPnIVMm25U3BZEEIYRIu0JDzeiljh0hc2bCIqP5eOkhlu27zICS6Xh/zAzUb789/hxPMyVHGiMJQgiRdi1fDoGBxPTowYp9voxbe5LLt+/yQfMyvN24JNQ+bNPJ71IbSRBCiLTj3sI898yZw+mKNRl0yIUjVw5QoUAWxnasRN2Slqktype3fRlTEUkQQoi0Yfp0ag0fblZ0K10aLl3i1L4TdO49CQLDmdSpCi9Xyo+Tk+ONRkoqSRBCiLRhxgzc/f3hxRfBx4eTsxfzv06jcMqQngVv1qJEbsfsaH4WMsxVCJH6+fvDzp1cb9AArl7l1Ks96HyzAE6urizoV1eSQxJJDUIIkfpZpu++0K0bdwcMocvmWzhFRbKwEhSX5JBkUoMQQqR+K1dCoUL4FixO13OZCM+Snfmnl1P8fwms0iYSTWoQQojULSwM1q3jdvc+jNsdxs0IJ359qx5lCr9s75KlelKDEEKkbt7eXHHOwKt5X+BqiGb6615ULZzd3qVKE6QGIYRI1Y6t9qbH6+MJjXbmPS/X2HscxDOTGoQQImXasQMOHnzsLjvPBvCqqxcqnRu/96tDuZxJXNJTJEgShBAiZercGVq1Mn0MCbgQEMIbc3aSJ+gGS0uF8lzeLDYuYNonCUIIkfL4+cG5c2aa7h9/fOjt4PAo3pi3GxURwezfPyd/25fsUMi0zy4JQin1rlLqiFLqsFJqgVLKXSlVTCm1Qyl1Sim1SCnlZo+yCSFSAB8f81ikCIwaBXfu3H8rJkbz7qL9nLkewpQtP1O4UmnIm9dOBU3bbJ4glFIFgIGAl9a6AuAMdALGABO01qWAW0BvW5dNCJHMDhwwK7o9LR8fcHMzU3XfvAnffgtAVHQMI1ceYf1Rfz4t5UJdnzXQt28yF1rcY68mJhcgvVLKBcgA+AFNgCWW9+cCbe1UNiFEctizB6pUge+/f/pjt20DT0+oVQteew0mTOD2uUv0nLOLuT4X6F2vGD3+nmmWEO3QIfnLLgA7DHPVWl9WSo0DLgJ3gXXAHuC21jrKspsvUCCh45VSfYG+AB4eHnh7eyepHMHBwUk+NjVzxLgdMWawf9yF58+nOBD52WdsL1mS6DirsrkEBhKVOXOC6zyryEjq79rF5XbtOOPtjXurl3H57xADf9jKtXRZ6VnBjSaRJ9BLl+Lbti1nduy4f6y9Y7YXq8WttbbpD5Ad2AjkBlyB5UA34HScfQoBh550Lk9PT51UmzZtSvKxqZkjxu2IMWudAuJu2lTr3Lm1Bq2HD4/dvnWr1unSad2zp9YxMQ8ft3271qCvL1iiZ2w9q5tP2KyLDF2la/Sfq/dOnGX2GTPGnPfo0XiH2j1mO3nauIHdOhHf1/ZoYnoeOKe1vq61jgSWAnWAbJYmJ4CCwBU7lE0IkRzCwuC//6BLF3j1VRg/3sy4eu4ctGsHLi4wezZMmxbvML87d5n3z1E6dxpFjQPufLnqKOlcnPjy5bKsv7iMqu+/AWvXwvTpUL8+lC1rpwAdgz3upL4I1FJKZcA0MTUFdgObgI7AQqA7sMIOZRNCJIdt20ySaNoUypSBP/6AoUNh926IijKP772HfucdNniUZUNMNrafvcm5GyFAHkpmi+btxiVpVSk/ZfJmNuf8bR7UqwcvvwyRkfDZZ3YN0RHYow9ih1JqCbAXiAL2AdOA1cBCpdRXlm0zbV02IUQy2bgRnJ2hQQPIkgX69IGffzbb1q6F557j6tSZDBs+lw3bg8icLpSaxXPSpWZhGr7ViVKVSsALfeKfM3NmWLUKataE8HDo2NE+sTkQu8zFpLX+DHgw/Z8FatihOEKI5LZhA1SvbpIDwIgRsHUrvP8+ofUb8ofPecb+fYLIQpUYvnEW3Z39cfl4PVy/Dkd3wZtdEz5voUKwaxcEBoK7u83CcVQyWZ8QInkFBpov8Y8+it2WPz87V21l4a6L/P3VP4RGRFOreA7GdKhEkSohpp/i9dehfXuzf506jz5/gQLmR1idJAghRPLavBmio03/g8WK/ZcZtHA/mdO50LpyftpVLUCNYjlQSsErr8C4cTBkiDk2fXqoXNmOAYh7JEEIIZLXxo2m+ad2bQD2XLjJB78fpGaxHMzpWYP0bgnMuPree2aE05QpZnSSq6uNCy0SIglCCPFstAZfX8iXzwxf3bAB6tYFd3cu3Qyl77w95M/mzk9dPRNODmBumJs0yXRi169v2/KLR5LZXIUQSXfjBrRtC4ULm2kvmjSBQ4egaVOCwiLpNWcXkdExzOxRnewZnzD/prOzSRIyOinFkBqEECJpNm6Ebt1MkvjkEzPjqo8PZMtGTOvWvL/4AGdvhPBLrxqUyJ3pyecTKY4kCCHE01u5Etq0gdKlzb0JVavGe3vKhlOsO3qez14uRx1ZAjTVkgQhhHg6WsMXX0CpUmbG1owZ47298bg/4/85SfuqBehRp6h9yiiShSQIIcTT2brVTJXx00/3k8PNkAh2nb/JznM3WbzrEuXyZeHr9hXNMFaRakmCEEI8ne++g1y5zI1twMx/zzFq9VFiNKRzcaJ60Rx806Ei7q6PGLEkUg1JEEKIxDt50vQ/DB+Odnfn+39OMeGfkzQr58GbDYpTsWBW0rlIYkgrJEEIIRJvwgRwc0P378/oNceZtuUsHaoVZEyHirg4y6j5tEYShBAivuBgCAkBD4/422/cgDlzCO7Wk6EbLrP6oB+v1y7C5y+Xx8lJ+hrSIkn5Qoj4Bg6ESpXMpHtxTZrEyUx5aF2sHWsO+TH0xecY2VqSQ1omNQghRCytzXoN164RM34CB3oNZN/F21y+fIMrx8C75/dkjHHi1z61qF0ip71LK6xMEoQQItaZM+xRWVjY+n9sulmKGz9uAyC9jiJ/9gI8XzI7wzvVIE8WWYvBEUiCEEIA4HsrlG8W7mdV13FkcdY0OrSFpqVzUadbK3J5VkK9+Sb0fcvexRQ2JAlCCAfleysUnzMBXLoZyrmAUNYduYqKcGXQgRW8+ecUMrz1B0z+HvatNdN3Dx9u7yILG5MEIYQD8j5xjQG/7SM4PAonBfmzpadVpXy8P6I7+SuVgXSu8Pnn8NtvZlK+ESMgb157F1vYmCQIIRzMLz7n+ezPIzyXAb6r7EKJ1s1wc3GCs2fhxEF4p6/ZsUgRGDoUfv0V3n/frmUW9iEJQggHEROjGb3mGNO3nqNpiWx8P7QNGXW0SQzZs5vlPgEaNYo96IsvTO3BRb4qHJHcByFEWhYRAV27Ej1zFh8tPcj0refoXrsI086uJuPtm2YNh2++Mft6e5s5lsqVi38OSQ4OSz55IdKyYcOIXLCQd4MLsuo5DwY2LcW7VbKjuv4Ar71m1n7+/nt45x1Tg2jY0Cz/KQRSgxDCvm7fNlNnW0GOHTsImzCJfgN+ZNVz9fnk6Greq18YNX68mUpj+HDThBQTA336wIULJkEIYSEJQgh7GjUK6tY18x8lpytXKPLtd7zRfQz/pC/Al6UUfVdOhbfegh9+gE6dTFNS0aLQv7+5exri9z8IhycJQgh78vEx/QR79iTfOWNiCHm9J32bv8+/uUsztmMluvVuAW+/DXPnQmho/Hsahg2DzJkhRw4oXz75yiFSPemDEMJeoqJg3z7zfPv2ZGveuTljDn3yNWN//jJM7FSFNlUKmDfGjoUdO6BGDShbNvaAXLlg9mwICgIn+ZtRxJIEIYS9HDtm/poHkyCSwblTvvTcp7mSrxT9q7jHJgeADBlMgkioE7pDh2S5vkhb7PLnglIqm1JqiVLquFLqmFKqtlIqh1JqvVLqlOUxuz3KJoTN3OucrlHDJAitn+l0O8/dpN30nQS6pmfBSwXxyuf68E5OTjJKSSSaveqTk4C/tdbPAZWBY8BHwAatdSlgg+W1EGnXrl2QJQt06wZXr8LFi099iugYzT9H/ek5eyev/exDjlvXWJbhFJ6NPa1QYOFobN7EpJTKAjQAegBorSOACKVUG6CRZbe5gDcw1NblE8Jmdu0CT0+oU8e83r7dTG8BcP06/PKLWbznETeq+ZwJYMjvB7h8+y55Mrsx4OQ/9Dn0N1kPJGOHt3Bo9uiDKA5cB2YrpSoDe4BBgIfW2g9Aa+2nlMqT0MFKqb5AXwAPDw+8vb2TVIjg4OAkH5uaOWLcKTFmFRFB/f378e3YkXM3b1IvXTquLFnCGcsynyUnT6bgH39w/OpVrrZoEe/YGK1ZfTaSpaci8cioeKdqOjos/omiy/7g4OjR3Ny5E0iZcVubI8YMVoxba23TH8ALiAJqWl5PAr4Ebj+w360nncvT01Mn1aZNm5J8bGrmiHGnyJh37dIatF682LyuX1/rWrXM86AgrbNkMe8XKaJ1ePj9wy4GhOjus3boIkNX6QG/7dVBYZFar1hh9h00KN4lUmTcVuaIMWv99HEDu3Uivq/t0QfhC/hqrXdYXi8BqgH+Sql8AJbHa3YomxC2sWuXeaxe3TzWqgV790J4uGlaCgw0021fuAAzZ3I7NIKvVh2l6Xeb8TkTwJdtK/B9pypkun4VevaEqlVhzBi7hSPSJps3MWmtryqlLimlymitTwBNgaOWn+7AN5bHFbYumxDJytsbIiOhWbOH39u929x/cK/PoVYt+PZbc1/E5Mmmb2LECPjnH9bPWcn7/sUJCo+iY7WCvNusNPmzpYfoaOjSxSSVhQshXTqbhifSvkQlCKVUaWAqpp+gglKqEtBaa/1VEq/7DvCrUsoNOAv0xIyoWqyU6g1cBF5J4rmFSBn694dr18DX16zIFteuXeDlFTvktFYt8/j113D0KMyZgwZm9R7BV8fCqBgVzNhBzXgub5bYc/z1F2zZAtOnQ+nSNglJOJbENjFNBz4GIgG01geBTkm9qNZ6v9baS2tdSWvdVmt9S2sdoLVuqrUuZXm8mdTzC2F316+bG+ECAuD33+O/FxoKR47ENi8B5M8PhQrBypWQKxdRHV/hsz+P8OXxCJrfPsOiWYN5LrNz/PNMm2ZWeeve3frxCIeU2ASRQWu984FtUcldGCHSjH//NY8ZM8LUqfHf27fPzKAaN0HA/VpEWJ83eWvJUeb5XODNBsX5sXsN0vv5wsSJsfteumRqEL16mSm7hbCCxCaIG0qpEoAGUEp1BPysViohUrstW0yz0ogRZkK+e3MuQWwHtZcXUdExhIRb/tZq0oQ7WXLwet6mbDjuz8jW5fm4RVmc6tSBNm3MXEo3bph9Z8wwd16/8YZt4xIOJbGd1G8D04DnlFKXgXNAV6uVSojUbutWqFkT+vY1o5GmTjVNQufOse2XP1n0v5GcXHyaM9f3ExEVQ4ncGalcsCZHh//OmWt3+b5TVV6unD/2fF9/DRUrmsexY02CaN7cTNcthJUkKkForc8CzyulMgJOWusg6xZLiFQsMNDUGIYNg2zZzEijX38levC7TBo6hR+ef5ec6V2okCUd9UvlIqObC4cu32bLqRtERmtm9ahO/VK545+zXDnT1zBlChQrBleumOdCWFFiRzF9DYzVWt+2vM4OvK+1/tSahRMiVfLxMX0M9eub1/36EfDrYgZ9tZx/y79Eh8Lp+KpPY9K7xe901loTo8HZ6RGT6Y0cCb/9BoMGQb580LKllQMRji6xfRAv3UsOAFrrW0CLx+wvhOPasgWcnaF2bQACy1Wkc++J7MpfhjHlXBjXr+lDyQFAKfXo5ABmlNOAAabvoXdv6ZwWVpfYPghnpVQ6rXU4gFIqPSB35QiRkK1bzY1umTIRGR3D27/u5WzmPMxt4kHd572e7dzDhpkb4wYOTJ6yCvEYiU0Q84ENSqnZmJFMvTAzrgoh4goLM4vyvPMOWmtGrDjC1lM3GNuhEnWrF3r282fPbtaUFsIGEttJPVYpdQgzLYYCvtRar7VqyYRIDWJizFDTXLngs8/M2tIREdCgAdO3nmXBzov0b1SCV5MjOQhhY4mei0lrvQZYY8WyCGF/oaHwzTdmxFCJEk/ef8IEmDXLPF+2DKpVA2BN7rJ8veI4LSvmY8gLZaxYYCGsJ1Gd1Eqp9palQO8opQKVUkFKqUBrF04Im5syBb780oxAOnLk8fseOgSffALt2sE//5iaw6JF7Kv3EoNXn6Zq4Wx892plnB7X8SxECpbYGsRY4GWt9TFrFkYIu7p7F777zkyid/kyNGwIa9eaDucHhYcT1KM33lVfYGObD3G9lZ4qc9dR5O9lDHIuT54s6Zj+uhfurg+PVhIitUhsgvCX5CDSvFmzwN/fTJ1dqBA0bQpNmsCCBWBZ1S0sMpp1R66y/Ld/+LfxJ0S4uJLjYhAxOpDFoZFABbK4uLCwRw1yZZKBfiJ1S2yC2K2UWgQsB8LvbdRaL7VKqYSwtYgIM4VF3bqm5qCUmXCvZUto2ZILAz/kp7qdWHXIj6BITb5gzevpLtJ8QFeqFc6Ok4ILAaEcvHyHCvmzUDx3JntHJMQzS2yCyAKEAi/E2aYBSRAibZg/Hy5eRE+dypErgUREx5DFPRsxKzcw84elLFF5cdl5nlbH/6XD+R3UGtANpzfeNjfEWRTNlZGiuTLaMQghkldih7n2tHZBhLCb6GgzcqlaNX7OUJpvfvg33ttu6QrRLS/0/3kkeWpUgdnLIHfuR5xMiLQjsXMxuQO9gfLA/aWxtNa9rFQuIWxn7Vo4dYoVPy3lm79P0LJSPjp6FiTwbiShEdE0KpObfFnTw4CX7F1SIWwqsU1MvwDHgebAF0AXQDqtRdqwYgXbnqvFkIvpqFU8O+NfrUw6Fxl9JERiJ+srqbUeDoRorecCLYGK1iuWEDYSE8Ph//bzZqsPKZYrIz9385LkIIRFYhNEpOXxtlKqApAVKGqVEglhQ8c27qBr08FkcXdmTs8aZE0vM6QKcU9im5imWdaA+BT4E8gEDLdaqYSwgZP+QXT5xx/3qAgWdK9P/mzp7V0kIVKUxCaIDZY1ILYAxQGUUsWsViohnsWaNWbqi3HjzP0MQEyMZsqm08zfcYHM7q7kzOjGqWvBuITdZcHJPyhcsrudCy1EypPYJqY/Eti2JDkLIkSyiIiAt96C8eNh2zYAgiI0Pebs4rv1JyntkZnSHpnQGkpldeW3Xz6kWLN6di60ECnTY2sQSqnnMENbsyql2sd5KwtxhrsKYTMXL4KLC+TPn/D7s2axNl1+5nTuR4YFB8l8KQNbjt8lOCqMr9tVpHONQihLrYIff4SbvtC6te3KL0Qq8qQmpjJAKyAb8HKc7UHAG9YqlBCP1K4dZMkCmzY9/F54OCcnz2JQm0/I6RxNVv/LnDpznSxuinlv1KFCgazx91+5EkqWhDIyHbcQCXlsgtBarwBWKKVqa619bFQmIRIWGAj79pnpLYKDIVP8+Y7uTp/FgJo9yOTuyrLXq5KnfGfo1w/vdu1Mcvj2WzOVd8OGZvK9jRvNGs9KpuMWIiGJ7YNop5TKopRyVUptUErdUEp1tWrJhHjQzp2gNURFmYn04goL44tN5zmZuwjju1YnT6ki0KkTzJyJS3AwTJwIH34IFSqYdRz69zf9FdK8JMQjJTZBvKC1DsQ0N/kCpYEPrFYqIRLi42P+2nd1hQ0b4r31x6SFLCjVgH7FXGhQJo/Z+O67EBJChU8/Nc87dIAtW+DcOZMkli6FBg3sEIgQqUNiE8S9u4daAAu01jef9cJKKWel1D6l1CrL62JKqR2WlesWKaXcnvUaIo3Ztg3KlYM6dUzzkMXMrWcZcjMntW6f573ez8fuX7UqNG5MtgMHoFUr+O0308GtlKlJtGsnzUtCPEZiE8RKpdRxwAvYoJTKDYSx/SlDAAAgAElEQVQ947UHEX8+pzHABK11KeAWZnJAIYyYGNi+neg6dbjbuCl63z6ibwQwcuURvlx9jOYnfZhTLR2uD06TMWkSF/73P/j9d3CTvzmEeBqJShBa64+A2oCX1joSCAHaJPWiSqmCmPmcZlheK6AJsfdWzAXaJvX8Iu0JOnSMn0s3oU6eVpS9W40SQ5ZTYbwPs/87T6/ws0xZ9z3unV59+MCKFTn3xhvgLqOyhXhaSmv96DeVaqK13vjAPRD3JXVFOaXUEmA0kBkYAvQAtmutS1reLwSs0VpXSODYvkBfAA8PD8+FCxcmpQgEBweTKZPjrfqVkuJ2Cg+n0gcfcMvLiwtdu4LTw3+vRMdoVp2NZN2pu4QoFypkiOS5fO54/LaI66XKkL1eFT7o35Ebdepw/JNPErxOSorZlhwxbkeMGZ4+7saNG+/RWns9ab8n3QfRANiIuQdCA+qBx6dOEEqpVsA1rfUepVSje5sT2DXBzKW1ngZMA/Dy8tKNGjVKaLcn8vb2JqnHpmYpKu5du+DQIbIdOkSxK1fMqm65ct1/2z8wjHcW7GPnuVCaR1yj/18/U/nYTpNIln0PW7dD02IQEkLeoUPJ+4i4UlTMNuSIcTtizGC9uJ/UxBSklHoPOBzn5whwyPI8KeoCrZVS54GFmKaliUA2pdS9hFUQuJLE84vU4rDln9DHH4O3N1Srdn/b1lPXaTFpK4d87zDhtcr8vG4ilUt6xNYymjSB48dhzBgoVAgaN7ZPDEKkYU9KEJkwzUCeQD8gH5AfeAsol5QLaq0/1loX1FoXBToBG7XWXYBNQEfLbt2BFUk5v0hFDh82fQNffgn//Qfh4UT2f5sxa47x+qyd5Mjoxp8D6tKuWCY4ehRq1449tmlT87hvH3TvnmDzlBDi2TzpTuqRAEqpdUA1rXWQ5fXnwO/JXJahwEKl1FfAPmBmMp9fpDSHD5thq87O4OnJpWFfMnBXEPs2n6VT9UKMeLkcGdxczJKgED9BVK4MOXLAzZvQo4ddii9EWpfY6b4LAxFxXkeQDAsGaa29AW/L87NAjWc9p0hFDh9mTavubP7jIPsv3eakfwEy5r7LD4d+5+XRs2PvUdi2zdQQasT55+HkBJ07w+XLUKKEfcovRBr3NGtS71RKLcN0HrfDDEUVImlu3WJdhkL0y16XrIevUrlQNl4on5dXzvlQaMJcWN7G3MgWEWHWdqhYETJnjn+OyZPtU3YhHESiEoTWepRSag1Q37Kpp9Z6n/WKJdK6kAOH+bzZm5TJAKuGPY+rs6UPIao4jB8Nw4dDnjxmbYfDh01ntBDCphJbg0BrvRfYa8WyCAfy/TZfrmTJw5LmxWOTA5ipMEaONBPt1atnRij9+Se8/PKjTyaEsIpEJwghksvxq4HMuJOJ145txGt0i4d3eOUVWL4c8uUzyeLBpiUhhE1IghA2dTcimmHLDpMlKoyPbu9PeLI8JydYsMD2hRNCxCMJQthEVHQMi3f7MmnDSfwDwxn/3y9krySjj4RIySRBCKuKio5h1UE/vt9wirM3QvAskp0fmhehxpiV8L9J9i6eEOIxJEEIq4iJ0Szbd5nJm05z7kYIz+XNzPTXvXi+bB7UvcV+ype3byGFEI8lCUJYxfj1J5m86TRl82Xhp67VeKFcXpycLP0N9+ZgqvDQZL1CiBREEoRIdmsO+TF502le8yrENx0qoh7siD5yxMzamiePfQoohEgUmeFMJKsTV4N4//cDVM3jzhd39qDu3Hl4p8OHTe1BlvsUIkWTBCGSzbWgMPpO/5eMoUH89Hkn0vXpBWXKwNy5ZslQAK1jE4QQIkWTBCGeWXSMZu628zQdvR6/O2H8tGosHm/3gfXroXhxM9tq9erQqJF5HRwsCUKIVED6IMQzOX41kCG/H+Dw5UDqXT3ByLPrKbFrY+zdz02awJw5MGkSREdDnTrQrRu8msD60UKIFEUShEiyFfsvM/SPg2R2d2Vy7hu0HPMBas2a+FNjODlBr17mRwiRqkiCEI91524kfx/2Y8X+K1wNDKNG0RzULpGTA5fuMOu/c9QomoPJHcqRp2oPUzto3tzeRRZCJBNJEOKRfthwih82nSYiKoYiOTNQPFdGVh/yY+GuSwD0rFuUT1qUxfXHKWbhnnnzZGSSEGmIJAiRoNUH/fhu/UleLJ+XtxqVoHLBrCiliI7RHLlyh6gYTbXC2SEkBL7+2nRAN2li72ILIZKRJAgBvXuDvz989BHUq8fpa0F8sOQA1Qpn4/vOVXFziR3s5uykqHTgPxg61BwTEGCGri5ebMcAhBDWIAnC0R06BLNmgasrrF5NcKPnefP598ng5syPXTzjJYf7xo+HGzegY0dzN7SXF9Sv//B+QohUTRKEo/vhB3B3h5Mn2btwNaOOhXP+TgTze9cgb1b3h/cPCIDNm00NYtQo25dXCGEzcqOcI7t5E+bP5+jr/ei94SrtAwpxrkBJxq2eQO1taxI+ZuVKcz9D+/a2LasQwuYkQTiwuzNmMbpmJ17O0ZTdF27xQfMybB32Au3cbsOYMbHTY8S1dCkULgzVqtm+wEIIm5IEkZbdugXPPw/DhpnRRha3QyP4++BlXryQi59rduAVr0Js+aAxbzcuSUZ3V9N8dOIErFgR/3zBwbBuHbRrJ8NZhXAA0geRVmkNffrApk2wYQMRv8xn9MCJ/BWaDv+/1wNQJCqS38pEUKdDpfjHduxokso330DbtrHJYM0aCA83CUIIkeZJgkirfv7ZNAd9+y0RNWvRf/Z2/rnhxvOnt9M9xI8qV05S7c4l3KecfPhYFxcYMgT69zcd0o0ame3LlkHu3FCvnk1DEULYhzQxpUWHD8O770Lz5kQMHEz/s+78k6csXxYMY3i6M/TPFUadnM64f/O1SQYJ6dHDDGEdPBgOHjQ1h1WroHVrcHa2aThCCPuQGkRaEx7O9W69Of5cbc4MHsu6ubvZdiaAL9uUp1vtonhXyEmRezWCx0mfHqZONc1UVaqYWkRQkIxeEsKB2DxBKKUKAfOAvEAMME1rPUkplQNYBBQFzgOvaq1v2bp8qZnWmpnjFjC62TCinZzB+xKZ3V34qm0FutYq8vQnbN8eGjc29zv88ANkzQpNmyZ/wYUQKZI9ahBRwPta671KqczAHqXUeqAHsEFr/Y1S6iPgI2CoHcqXKoVGRPHhwr2sCsrNC3dO0/PDrpTIk5HcmdI9vCb008ieHcaNg4EDzSimdOmSr9BCiBTN5glCa+0H+FmeBymljgEFgDZAI8tucwFvJEE8kdYa7xPXGb3mGKf9g/jQey79Jr6PKpEzeS9UuHDynk8IkeLZtQ9CKVUUqArsADwsyQOttZ9SKo8di5biaa1ZedCPHzed5vjVIPJndmPuytHUL50HatWyd/GEEGmA0lrb58JKZQI2A6O01kuVUre11tnivH9La509geP6An0BPDw8PBcuXJik6wcHB5MpU6akFd5OMp04QdnRo9k/7jvmXM/E+gtR5M+oaFHclf+tnEWxxYvYNXMmocWKPfIcqTHuZ+WIMYNjxu2IMcPTx924ceM9WmuvJ+6otbb5D+AKrAXei7PtBJDP8jwfcOJJ5/H09NRJtWnTpiQfazcjR2oNetyIWbrI0FX68z8P6+joGK2vXNHa3V3r119/4ilSZdzPyBFj1tox43bEmLV++riB3ToR39U2vw9CmR7TmcAxrfX4OG/9CXS3PO8OrHjwWIe3bx8/1uzID+F56FS9ECNalcPJSZnhqOHhMHy4vUsohEhD7HGjXF2gG9BEKbXf8tMC+AZoppQ6BTSzvBZx/HUDxjbqQZsj3owqeNeMTgoPN3dNt2wJJUvau4hCiDTEHqOY/gUeNe5SBtk/wh2/63xWuT0VCWbcxqk4F7kLNWualdyuXTPDUIUQIhnJVBupxLilewnIkJWvvbLh2vplWLDA1B6+/x7KljWztgohRDKSBJEK7Lt4i/mXIum+dxUVG1SF7t3NYj/Dh8Pu3fDOOzL9thAi2UmCSOGiomP4ZNlhPKLv8v6pf8wEes8/D/nywbffmukvunWzdzGFEGmQJIgUTGvNV6uPccwvkM/3/0GmiuXMGy4u0KWLed6rFzjguG8hhPVJgkihtNaMWn2MOdvO06tmIZpvWgJVq8bu0L+/mWF18GC7lVEIkbbJdN8pkNaab/4+zox/z9G9dhGG5wtFRUebabfvKVbMrBYnhBBWIjWIFCYwLJIPlhzk581n6VqrMJ+3Lo/av9+8GbcGIYQQViY1iBRk88nrfPTHQfwDwxjQuCTvNSttbobbt890Rj9mjiUhhEhukiBsLDwqmqhoTcZ0sb/6O3cj+Xr1MRbtvkSJ3BlZ2r8uVQpliz1o/37TvCRDWYUQNiQJ4llERcGLL8KAAdC27WN31VqzfP9lRq0+RmBYFK0q5aPL3XPc8PVnOCUJCIngzYbFeff50rjv3A5Dp5j7HEqXNmtC9+1ro6CEEMKQBPEstm2DDRtMonhEgoiMjuHgvGWMverOjjtQuVA2mufPwor9l1ka7gYUomxeV2Z2r07FglnNQdOmmTul//gDevSA0ND4HdRCCGEDkiCexcqV5nHLFvDzMzevWSzadZE/D1xh34VbhEa6kzU8mNEdvXitdnGcnBQfZ7zGyg8nEq2ceO2nz3G9lxwAtm+Hhg0hf36TLEA6qIUQNiejmJ7FqlVmBlWtzV/7mKakcWtPMPSPQ1wPCueVGD8mr/iGLVN70/nEZjM9N5Dplzl0PrmFrvvX4LpzR+w5AwLg5EnTdPXbb7B6NXzyCZQvb48IhRAOTGoQSXX6NBw/DpMmwfTpsGgR+u23Gb3mONO2nKVT9UJ83bw4TkU7QJ064FYCJk6Enj0hMBCWLDHNRytXgo+PmU8JYIclWdxbNrRFC/MjRBoSGRmJr68vYWFhyXrerFmzcuzYsWQ9Z2rwqLjd3d0pWLAgrq6uSTqvJIikWrXKPLZqBXfucHP0OL6dv4MFRwLoVqsII1uXx2nKZFMjGDrU1Ap69TJ9FmfPQliYSRY3bpgmpXu2bwcnJ6he3T5xCWEDvr6+ZM6cmaJFi5qh3MkkKCiIzJkzJ9v5UouE4tZaExAQgK+vL8WSOERempiSatUqKFcOv5z5+KJIY+q+NYsFRwJ4s0FxvmhTHqfoKPjuO1N7qFsXOnc2E+1NnAizZ0O5ciYJ1KoF586Bv785r48PVKoEGTPaNz4hrCgsLIycOXMma3IQ8SmlyJkz5zPV0iRBJMWdOxw9eoH3Xn6f+mM2Mfd4IC9dP8r6PdP4uEVZ84/+99/hwgVTewBwdzfzJ61ebWoJvXqZ+xpq1zbv+/hAdLRpYrrXvCREGibJwfqe9XcsTUxxxcSYRXjSp3/kLjdDInhvyma8X59IBmfoVqsIvesVo+DU/TD/TzP0df16mDLFLOTTqlXswf36wejRZlhs165mW7Vq4OpqkkbJkhAUFJs0hBDCjiRBnD4NgwaZxwsXTJI4cgRKlXpoV//AMLrO2MHFO5oPdy2my58/kzWTu3nzlVdg2DDTnATQrJlJBk5xKml58pib34KCwMPDbHN3N0NYfXxi15SWGoQQaZrWGq01Tk4puxEnZZfOFpYvh7/+gooV4c03ITIy9v6GOC7dDOWVqf9x5UYQc1aPpX+BmNjkACahfPqpSQBnz8K6deDp+fD1hg2Db76Jv612bdi1C/79F3LkSDA5CSGSV0hICC1btqRy5cpUqFCBRYsWsWfPHho2bIinpyfNmzfHz88PgEaNGjF06FBq1KhB6dKl2bp1KwBHjhyhRo0aVKlShUqVKnHq1CkAxo8fT4UKFahQoQITJ04E4Pz585QtW5b+/ftTrVo1Ll26ZJ/An4LUIE6fhpw5zbBTMM1D69Zxp987zPr3HL637nIjMJRDJ64QFRnF/MUjqBpyFaZ9+/C5vvwyaWWoVcsMl/39d7PGg7TNCkcyeLCZbywZpI+OBmdnM/OA5Yv5Uf7++2/y58/P6tWrAbhz5w4vvfQSK1asIHfu3CxatIhhw4Yxa9YsAKKioti5cyd//fUXI0eO5J9//uGnn35i0KBBdOnShYiICKKjo9mzZw+zZ89mx44daK2pWbMmDRs2JHv27Jw4cYLZs2fz448/Jku81iYJ4vTp2KYdgBde4OqvS+g+9T9OXQ/BI4s7uSKC8Tyzj/dyhVB2zg/mS9zNLfnKcK/PITRUmpeEsJGKFSsyZMgQhg4dSqtWrciePTuHDx+mWbNmAERHR5MvzuwI7du3B8DT05Pz588DULt2bUaNGoWvry/t27enVKlS/Pvvv7Rr146MlpGI7du3Z+vWrbRu3ZoiRYpQKxX9H5cEcfo01Kt3/+Wpus3oHlGZwIAQ5vWqSb1SuUxH84EDpo/CGm2GhQtD3rxw9ap0UAvH84S/9J/G3ae4D6J06dLs2bOHv/76i48//phmzZpRvnx5fHx8Etw/Xbp0ADg7OxMVFQXA//73P2rWrMnq1atp3rw5M2bMQGv9yGtmTGXD1x27DyI8HC5evF+D2H42gI5HXIh0cmFR1B6THAICYO1a6NTJOskBYoe7KiU3yAlhI1euXCFDhgx07dqVIUOGsGPHDq5fv34/QURGRnLkyJHHnuPs2bMUL16cgQMH0rp1aw4ePEiDBg1Yvnw5oaGhhISEsGzZMurXr2+LkJKdY9cgzp0z8yiVLMnvuy/xybJDFM6RgTnb/6CQ/wUYM8LMsRQVZW50s6YhQ8wIqKxZn7yvEOKZHTp0iA8++AAnJydcXV2ZOnUqLi4uDBw4kDt37hAVFcXgwYMp/5h50BYtWsT8+fNxdXUlb968jBgxghw5ctCjRw9q1KgBQJ8+fahater9ZqlU5d5wq9T44+npqZNq06ZNWq9cqaNReuzMDbrI0FW6y/Tt+nZohNZff601aH31qtaNG2tdurTWMTFJvlZKsmnTJnsXweYcMWatU3bcR48etcp5AwMDrXLelO5xcSf0uwZ260R8xzp0E9OVE+fp9tqXTDl5l841CjO7Z3WypneFF14wO8ybB97epvYgI4uEEA7GIZuYtNb8dzmSd64VJLpAPka3q0CnGoVjb0uvWtUMfR050jRBdepk3wILIYQdpKgahFLqRaXUCaXUaaXUR9a6zuSNp5l+KIIyIdf5e+fPdK5ZJP6cJU5O5k7okBAznvq556xVFCGESLFSTIJQSjkDU4CXgHJAZ6VUOWtcq121ArxWxo1Ff4+lcMGcCe90r5nJ2p3TQgiRQqWYBAHUAE5rrc9qrSOAhUAba1yoYPYMtCikcD5/Lv5NcnF16AADB0Lv3tYoghBCpHgpqQ+iABB3chJfoKa1Lubu72+m135UgsiSxUx/IYQQDiolJYiEhgk9dEuiUqov0BfAw8MDb2/vJF0s/enTAOwNDCQwiedIjYKDg5P8O0utHDFmSNlxZ82alaCgoGQ/b3R0tFXO++uvv7J3716+++67pz72rbfe4sUXX6Rt27aPPX+TJk3uT+0xYMAABgwYwHOJ7P98XNxhYWFJ/neQkhKEL1AozuuCwJUHd9JaTwOmAXh5eelGjRol6WKnli0DoNqrr8ZOve0AvL29ServLLVyxJghZcd97NgxqywNaq0lR93d3XFzc0vSuV1dXUmfPv1jj124cCFeXl6ULl0agLlz5z7VNR4Xt7u7O1WrVn2q892TkvogdgGllFLFlFJuQCfgT2tdLP3ly5Apk1mjQQjhkObNm0elSpWoXLky3bp1Y+XKldSsWZOqVavy/PPP439vKeA4/P39adeuHZUrV6Zy5cps27aN8+fPU6FChfv7jBs3js8///yhY7/44guqV69OhQoV6Nu3L1prlixZwu7du+nSpQtVqlTh7t27NGrUiN27dwOwYMECKlasSIUKFRh6b4VKIFOmTAwbNozKlSvTpEmTBMv6rFJMDUJrHaWUGgCsBZyBWVrrx0+E8gzSX75s+h/kBjgh7GrkyiMcvRKYLOeKjo7G2dmZcvmz8NnLj54iA8xaDqNGjeK///4jV65c3Lx5E6UU27dvRynFjBkzGDt27EPNSgMHDqRhw4YsW7aM6OhogoODuXXrVqLKN2DAAEaMGAFAt27dWLVqFR07dmTy5MmMGzcOLy+vePtfuXKFoUOHsmfPHrJnz84LL7zA8uXLadu2LSEhIdSqVYtRo0YxePBgpk+fzqeffvoUv60nS0k1CLTWf2mtS2utS2itR1nzWvcThBDCIW3cuJGOHTuSK1cuAHLkyIGvry/NmzenYsWKfPvttwlO1rdx40b69esHmJldsz7F/GmbNm2iZs2aVKxYkY0bNz5xMsBdu3bRqFEjcufOjYuLC126dGHLli0AuLm50cqypHGVKlWsMtdTiqlB2FR0NO5+fpIghEgBnvSX/tN4mj4IrXX8G2SBd955h/fee4/WrVvj7e2dYDNRQlxcXIiJibn/Oiws7KF9wsLC6N+/P7t376ZQoUJ8/vnnCe73YBkfxdXV9X75405BnpxSVA3CZi5dwikqShKEEA6sadOmLF68mICAAABu3rzJnTt3KFCgAPDojuKmTZsydepUwDRpBQYG4uHhwbVr1wgICCA8PJxVq1Y9dNy9ZJArVy6Cg4NZcm8VSyBz5swJjkKqWbMmmzdv5saNG0RHR7NgwQIaNmz4bIE/BcesQVjWjZUEIYTjKl++PMOGDaNhw4Y4OztTtWpVPv/8c1555RUKFChArVq1OHfu3EPHTZo0ib59+zJz5kycnZ2ZOnUqtWvXZsSIEdSsWZNixYolODw1W7ZsvPHGG1SsWJGiRYtSPc7aLz169OCtt94iffr08RYsypcvH6NHj6Zx48ZorWnRogVt2ljl/uEEqcdVYVI6Ly8vfa+n/6lMnQr9+4OvL1j+WnAUKXnoo7U4YsyQsuM+duwYZcuWTfbzWmuYa0r3uLgT+l0rpfZorb0SPCAOx2xiyp+f6/XqQZz1ZoUQQsTnmE1MbdpwJGtWGllrCVEhhEgD5BtSCCFEgiRBCCHsIjX3f6YWz/o7lgQhhLA5d3d3AgICJElYkdaagIAA3N3dk3wOx+yDEELYVcGCBfH19eX69evJet6wsLBn+kJMrR4Vt7u7OwULFkzyeSVBCCFsztXVlWLFiiX7eb29vZM8c2lqZq24pYlJCCFEgiRBCCGESJAkCCGEEAlK1VNtKKWuAxeSeHgu4EYyFie1cMS4HTFmcMy4HTFmePq4i2itcz9pp1SdIJ6FUmp3YuYiSWscMW5HjBkcM25HjBmsF7c0MQkhhEiQJAghhBAJcuQEMc3eBbATR4zbEWMGx4zbEWMGK8XtsH0QQgghHs+RaxBCCCEewyEThFLqRaXUCaXUaaXUR/YujzUopQoppTYppY4ppY4opQZZtudQSq1XSp2yPGa3d1mTm1LKWSm1Tym1yvK6mFJqhyXmRUopN3uXMbkppbIppZYopY5bPvPaDvJZv2v5931YKbVAKeWe1j5vpdQspdQ1pdThONsS/GyV8b3lu+2gUqras1zb4RKEUsoZmAK8BJQDOiulytm3VFYRBbyvtS4L1ALetsT5EbBBa10K2GB5ndYMAo7FeT0GmGCJ+RbQ2y6lsq5JwN9a6+eAypj40/RnrZQqAAwEvLTWFQBnoBNp7/OeA7z4wLZHfbYvAaUsP32Bqc9yYYdLEEAN4LTW+qzWOgJYCNhuFXAb0Vr7aa33Wp4HYb4wCmBinWvZbS7Q1j4ltA6lVEGgJTDD8loBTYAlll3SYsxZgAbATACtdYTW+jZp/LO2cAHSK6VcgAyAH2ns89ZabwFuPrD5UZ9tG2CeNrYD2ZRSSV5b2RETRAHgUpzXvpZtaZZSqihQFdgBeGit/cAkESCP/UpmFROBD4EYy+ucwG2tdZTldVr8vIsD14HZlqa1GUqpjKTxz1prfRkYB1zEJIY7wB7S/ucNj/5sk/X7zREThEpgW5odyqWUygT8AQzWWgfauzzWpJRqBVzTWu+JuzmBXdPa5+0CVAOmaq2rAiGkseakhFja3dsAxYD8QEZME8uD0trn/TjJ+u/dEROEL1AozuuCwBU7lcWqlFKumOTwq9Z6qWWz/70qp+Xxmr3KZwV1gdZKqfOYpsMmmBpFNksTBKTNz9sX8NVa77C8XoJJGGn5swZ4Hjintb6utY4ElgJ1SPufNzz6s03W7zdHTBC7gFKWkQ5umE6tP+1cpmRnaXufCRzTWo+P89afQHfL8+7ACluXzVq01h9rrQtqrYtiPteNWusuwCago2W3NBUzgNb6KnBJKVXGsqkpcJQ0/FlbXARqKaUyWP6934s7TX/eFo/6bP8EXreMZqoF3LnXFJUUDnmjnFKqBeYvS2dgltZ6lJ2LlOyUUvWArcAhYtvjP8H0QywGCmP+g72itX6wAyzVU0o1AoZorVsppYpjahQ5gH1AV611uD3Ll9yUUlUwHfNuwFmgJ+YPwDT9WSulRgKvYUbt7QP6YNrc08znrZRaADTCzNjqD3wGLCeBz9aSKCdjRj2FAj211ruTfG1HTBBCCCGezBGbmIQQQiSCJAghhBAJkgQhhBAiQZIghBBCJEgShBBCiARJghAikSwzpva3PM+vlFrypGOESM1kmKsQiWSZ02qVZeZQIdI8lyfvIoSw+AYooZTaD5wCymqtKyilemBm03QGKgDfYW5Y6waEAy0sNzGVwEw1nxtzE9MbWuvjtg9DiMSRJiYhEu8j4IzWugrwwQPvVQD+h5lOfhQQapk4zwd43bLPNOAdrbUnMAT40SalFiKJpAYhRPLYZFl3I0gpdQdYadl+CKhkmVW3DvC7mQ0BgHS2L6YQiScJQojkEXeun5g4r2Mw/8+cMOsUVLF1wYRIKmliEiLxgoDMSTnQshbHOaXUK3B/7eDKyVk4IZKbJAghEklrHQD8Z1k8/tsknKIL0FspdQA4Qhpc6lakLTLMVUAJhpgAAAA7SURBVAghRIKkBiGEECJBkiCEEEIkSBLE/9urAwEAAAAAQf7WC4xQEgGwBAHAEgQASxAALEEAsAQBwAquC+jhu9XuQQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107a46e48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sim_length = 100\n",
    "x, m = generate_simulation([0, 1], process_var=0.5, sensor_var=10.0, length=sim_length)\n",
    "\n",
    "plt.plot(range(sim_length), m[:, 0] ,'r', label=\"sensor\")\n",
    "plt.plot(range(sim_length), x[:, 0], label=\"calculation\")\n",
    "\n",
    "plt.ylabel(\"distance\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.legend(loc='lower right', shadow=False)\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of Kalman Filter\n",
    "\n",
    "To produce Kalman filter, we need to write two functions - predict and update You will see their code below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predict function\n",
    "\n",
    "As we already have seen above, prediction function consists of the follwing equations:\n",
    "\n",
    "$\\bar{\\mathbf x} = \\mathbf{Fx} + \\mathbf{Bu}$\n",
    "\n",
    "$\\bar{\\mathbf P} = \\mathbf{FPF}^\\mathsf T + \\mathbf Q$\n",
    "\n",
    "Prediction is used to compute prior - what we _predict_ will be the next state of the system. State of the system consists of mean of its variables and their covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(x, P, F, Q, B=0, u=0):\n",
    "    nx = np.dot(F, x) + np.dot(B, u)\n",
    "    nP = np.dot(F, P).dot(F.T) + Q\n",
    "    return nx, nP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Update Function\n",
    "\n",
    "Update function consists of the following equations:\n",
    "\n",
    "$\\mathbf y = \\mathbf z - \\mathbf{H\\bar x}$\n",
    "\n",
    "$\\mathbf K = \\mathbf{\\bar{P}H}^\\mathsf T (\\mathbf{H\\bar{P}H}^\\mathsf T + \\mathbf R)^{-1}$\n",
    "\n",
    "$ \\mathbf x = \\bar{\\mathbf x} + \\mathbf{Ky}$\n",
    "\n",
    "$\\mathbf P = (\\mathbf I - \\mathbf{KH})\\mathbf{\\bar{P}}$\n",
    "\n",
    "This function is used to choose middlepoint between our prediction and measurement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update(H, P, R, x, z):\n",
    "    y = z - np.dot(H, x)\n",
    "    A = inv(np.dot(H, P).dot(H.T) + R)\n",
    "    K = np.dot(P, H.T).dot(A)\n",
    "    nx = x + np.dot(K, y)\n",
    "    nP = P - np.dot(K, H).dot(P)\n",
    "    \n",
    "    return nx, nP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kalman Filter\n",
    "\n",
    "Using functions defined above, let's write Kalman Filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kalman(x, P, measures, R, Q, F, H, B=0, u=0):\n",
    "    xs, cov = [], []\n",
    "    for z in measures:\n",
    "        # predict\n",
    "        x = np.dot(F, x)\n",
    "        P = np.dot(F, P).dot(F.T) + Q\n",
    "\n",
    "        #update\n",
    "        S = np.dot(H, P).dot(H.T) + R\n",
    "        K = np.dot(P, H.T).dot(inv(S))\n",
    "        print(K.shape)\n",
    "        y = z - np.dot(H, x)\n",
    "        x += np.dot(K, y)\n",
    "\n",
    "        P = P - np.dot(K, H).dot(P)\n",
    "\n",
    "        xs.append(x)\n",
    "        cov.append(P)\n",
    "    return np.array(xs), np.array(cov)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters for our example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n",
      "(2, 2)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd8TecfwPHPkx2ZYsQKMWJHkBCqiFV7tVRVa7ZaqqW0tarVSVvVn9JWraKKGj+1ahNq116xBTFCEGSP+/z+uFd+aFREbm6S+32/Xnnde859zjnfbw73m/Occ56jtNYIIYQQD7OxdABCCCFyJikQQggh0iUFQgghRLqkQAghhEiXFAghhBDpkgIhhBAiXVIghBBCpEsKhBBCiHRJgRBCCJEuO0sH8DQKFiyofX19M7VsbGwsLi4uWRtQLmCNeVtjzmCdeVtjzvDkee/duzdKa13oce1ydYHw9fVlz549mVo2NDSUkJCQrA0oF7DGvK0xZ7DOvK0xZ3jyvJVS5zPSTrqYhBBCpMtsBUIpNUMpdU0pdeS+ed8opY4rpQ4ppZYopTzv+2y4Uuq0UuqEUqq5ueISQgiRMeY8gpgJtHho3jqgqta6GnASGA6glKoMvARUMS3zo1LK1oyxCSGEeAyzFQit9Rbg5kPz1mqtU0yTO4ESpvftgfla60St9TngNFDbXLEJIYR4PEueg+gNrDK9Lw5cvO+zCNM8IYQQFmKRq5iUUiOBFOC3e7PSaZbuk4yUUn2BvgDe3t6EhoZmKoaYmJhML5ubWWPe1pgzWGfe1pgzmC/vbC8QSqkeQBugif7/4+wiAJ/7mpUALqe3vNZ6CjAFICgoSGf2kja5HM56WGPOYJ15W2POYL68s7WLSSnVAhgKtNNax9330TLgJaWUo1KqNOAH7M7O2IQQIrf4dPOnrD+73uzbMedlrvOAHUAFpVSEUqoPMAlwA9YppQ4opSYDaK2PAguAY8Bq4C2tdaq5YhNCiNwqISWBTzZ/wtYLW82+LbN1MWmtu6Yze/q/tP8C+MJc8QghRF5w6sYpDNpAxYIVzb4tuZNaCCFykbCoMAAqFaxk9m1JgRBCiFwk7HoYCkX5AuXNvi0pEEIIkYscv3EcX09fnO2dzb4tKRBCCJGLhF0Po1Ih83cvgRQIIYTINQzawIkbJ6hYwPwnqEEKhBBC5Brno8+TkJIgRxBCCCEedDzqOEC2XOIKUiCEECLXyM5LXEEKhBBC5Bph18MomK8gBfIVyJbtSYEQQohc4viN49l29ABSIIQQItcIux6WbecfQAqEEELkClFxUdyIv4GfV0V+2HSa/RdumX2bUiCEECIXCLtuPEH927YUvllzgnXHIs2+TYs8UU4IIUTGXbubwKer1wFgl1qCGT2DaFzR2+zblQIhhBA5kNaa/4Yt5dylYszceosIfRh7eyc2D+6Mi6N9tsQgBUIIIXKgH3csY8C6jthqL1oU/YqCjnHcTqqYbcUBpEAIIUSOcjsumTGrwph0YAK2di54u7mz+tqbONo50r5C+2yNRU5SCyFEDqC1ZtnByzT9bjPz9x4n0X47vWq8wrG3DtCuQjvikuOoWrhqtsYkRxBCCGFh56JiGfXHEbaejsK/uActakfw+fZE3gx6HQ8nDxa/uJj1Z9dT16dutsYlBUIIISwk1aCZ9tdZvl13EkdbGz5tX4VuwaWoPW0gAd4B1CxaEwClFM3KNsv2+KRACCGEBVy4EceQhQfYfGEldh6r+KLVZ3T29+XA1QPsu7KP71t8j1LKojFKgRBCCDO5k3iHxccWc+z6McKiwriTeIcOFTtil1CfnzZdJNLmZ245rsEx2ZEu/23LkahRXI+7jqOtI92qdbN0+FIghBDCXD4J/YTxO8fjaOtI+QLlSU6FIWsHg7bFyd6DRMMthj87nA/qfcDA1QP5dMunAHSt2hUvZy8LRy8FQgghzEJrzeKwxbQo14J5HZfw/YazzNoRTmXny5Qvs5eY1NOMajiKBqUaADCz/UwalmrI6NDRDAweaNngTaRACCGEGey7so/zt8/TrvQ7NPtuK1ExiXStXZKhzZvhke+1f7RXStG7Rm961+htgWjTJwVCCCHMYPre+ShsWbKjCDVKODGtexABPp6WDuuJSIEQQogsFJ+UysSNp5i2dz758GdMh3p0rV0SWxvLXpGUGWa7k1opNUMpdU0pdeS+eV5KqXVKqVOm1/ym+Uop9b1S6rRS6pBSqqa54hJCCHPQWrPm6FWajt/MhM2hJKsIRjXpySt1SuXK4gDmHWpjJtDioXnDgA1aaz9gg2kaoCXgZ/rpC/xkxriEECJLnb8RS6+Zf/PGr3txc7KjbZ0LALwS0MnCkT0dsxUIrfUW4OZDs9sDs0zvZwEd7ps/WxvtBDyVUkXNFZsQQmSF+KRUvl17gmbjt7An/BYftq7E8refZffV1dQtUZfi7sUtHeJTye5zEN5a6ysAWusrSqnCpvnFgYv3tYswzbuSzfEJIcRjGbuTIvlsxTEuRcfTsHIydSvepmHZMly8E87+q/v5ptk3lg7zqeWUk9TpddDpdBsq1RdjNxTe3t6EhoZmaoMxMTGZXjY3s8a8rTFnsM68syPnq7EGfgtL4nBUKiVcFW/ViOfrs/2YfS4KAHtlfF5D0dtFs+33b668s7tARCqlipqOHooC10zzIwCf+9qVAC6ntwKt9RRgCkBQUJAOCQnJVCChoaFkdtnczBrztsacwTrzNkfOkTGRuDi44Gibj0kbT/PT9jM42NnwUZvKvBxcnBa/PUesIZalLy3lWuw1dkbsJL9Tfro9l31DZZhrX2d3gVgG9ADGml6X3jd/gFJqPhAM3L7XFSWEEJZU/5f6JKak4GczjtNXHGkXUIwPW1eisLsT7619j83nN/Nrx19pV6EdAK/V/OdNcLmV2QqEUmoeEAIUVEpFAB9jLAwLlFJ9gAtAZ1PzP4FWwGkgDuhlrriEECKjzt4M59TNUwBcZQDTOi7n1eAaJKcmM2n3JL7d8S1v1XqLV6q9YuFIzcNsBUJr3fURHzVJp60G3jJXLEII8aT2XbhF799/AaBO/qEcjf2JT3Z14XRcN6bum8qVmCs0LNWQ8c3HWzhS85FHjgohxH2u3I7ng0UHef7H7VxNOIiTbT62vPUZ615dy/W463y65VOqeVdjRdcVbOyxEQdbB0uHbDY55SomIYSwqOi4JH4KPcPM7eEYtKZvgzL8fj6cqi51sbe1J7hEMAfeOECqTqWcVzlLh5stpEAIIaxaQnIqs7aHM2nTaWISU+hYozjvNi2Pp0sqo746zIcVP0xrWzp/aQtGmv2kQAghrEaqIZWp+6bSrkI7iroWZdnBy3y9+gSXouNpVKEQQ1tWpGIRdwDWnVmHQRuoV7KehaO2HCkQQgirsfr0avqt7MeXf31FDedvORjuSOWi7nzdqRr1yhV8oO22i9uwUTbUKVHHQtFanpykFkJYjYXHFuNo68Kl21H8efUN3mrmzIq3nyW4jCcnb5wkOTU5re22i9vwL+yPu6O7BSO2LCkQQgirsOXUVeYeXIxdYi06lZyGp4sN3+x7iboz6uA2xo0Kkyrw+vLXAUgxpLAzYif1fKy3ewmki0kIkcfdik3isxXHmHtwFcmOdxjZsDsfN+3C8agAei/tjZOdE/2D+nMz4SYzD8ykednmVCxYkZikGKs+/wBSIIQQediGsEiGLj7M7fgkypY8wp0bzrzX4EUAKhasyPY+29PaphhSOBF1gjdXvknfmn0BrP4IQrqYhBB5zu24ZN5feIA+s/ZQ0NWBJf2fITwulBblWuDi4JLuMnY2dvz2/G8AjNsxjmJuxSjpUTI7w85xpEAIIfKMVINm9o5wqn/zJeOPNqRhtXMsHVCPGH2cy3cv80KlF/51+dL5S/NTa+MDLev51EOp3Pmo0KwiXUxCiBxp24VtRCdE07p86wy1/+vUdUZtj+dizHauu0xAG+4y98y7dDhVnB0RO7C3sc/Qul72f5m7iXcJLhH8tCnkelIghBA5zs34m3T4vQMJKQlcGXIFVwfXtM92X9rNwqMLeTv4bUp6lOT8jVg+XxnGumORFHJWlPH7ncjLMWztvZUha4fw4qIXcXd0p0mZJng6eWZo+28EvWGu1HIV6WISQuQ4IzaMICouipikGBYcXZA2X2vN26veZtyOcfhN9KPe5O40Gr+M7aejGNqiIi0q7mJLxDJGh4zmGZ9nWPPKGoKKBXEz/ibPV3zeghnlTlIghBA5yu5Lu5mydwqDggdRqWAlpu6bmvbZjogd7L60m7alB+Ca2pjtV38j3OEVUgsPYfut0Uw8M4GgYkF8UO8DANwd3Vnzyhp+bvMz3QO6WyqlXEu6mIQQOUaqIZV+K/tRxLUInzT6BB8PH4asHcKRa0eoXLAK7636AjvcOHAshBblutH1mc/YHbmM/Vf3s/HcRpIMScxsPxM7m/9/tbk7utM3sK8Fs8q9pEAIIXKMn/f+zL4r+5j3wjzcHd3pHtCd4RuGM3rDRK5dacyOO6sp7dSV6S/VJ6R8IZRStKkclLb8hk0bqFK4igUzyFukQAghcoTElEQ+3fwpDUs1pEuVLgDcjnWimGMDlpyYh7fdVWxtbNj85lh8PAunuw5bZZudIed5UiCEEDnC3MNziYyN5NeOv3LqWgwzt4ez4O+LYN8Yg816rqQu42X/l/HxLGHpUK2GnKQWQmQrrTXvrn6XuYfnPjBv/M7xlPGoxIwNrjz33RYW7Y3gxVo+7BwyiNKexgf1vFvnXUuFbZXkCEIIka1mHpjJf3b9B3sbeyoUqEBgsUCm7v6DI9eOUCBpEOGucXzQogIv1SqJl4vxec9jm45l96XdBBULeszaRVaSAiGEyDbXYq8xZO0Q6pSoQ8SdCF5c2IXW3rOZeuxz7Gy9GNPiDXrULY+D3YOdGy9WeZEXq7xooaitl3QxCSGyzaDVg4hNjmVGuxm8U30iZ6PPMTWsHwm2+xhWfyCv16/4j+IgLEf2hBAiW6w6tYp5R+bxVuD7TFgdxw9r7PFz6kWCzSGc7ZwZVOctS4coHiJdTEIIs7mdcJsN5zaw+vRqFh1bRGGnMizbVhs7dY33m1egz7M/8cbKZKoWqkqBfAUsHa54iBQIIYRZnLt1jqCpxnGQ8tm54ZhaDdv4l2kT4MMHLSpSzNMZgFkdZlk4UvEoUiCEEFnOoA30WdaHxJRkmhSYyKkIH/yLe/FJj6oElspv6fBEBkmBEEJkue93/sim8E0USnmbyCQ/PmtfgW7BpbC1se4H8OQ2FikQSql3gdcADRwGegFFgfmAF7APeFVrnWSJ+IQQ/3T21llGbhyJj7sPlQpWombRmgQUCXigjdaaGTt3M2Tt+zilVueVqr0Z3qoShdwcLRS1eBrZXiCUUsWBd4DKWut4pdQC4CWgFfCd1nq+Umoy0Af4KbvjE0Kk7/Mtn7Pw6EJsbWxJSjX+7fbtc98yuO5gAC5FxzPiv/tZEN4PG1sbFnaZSZsqAf+2SpHDWaqLyQ5wVkolA/mAK0Bj4GXT57OA0UiBECJHuBZ7jd8O/8YbgW8woeUEzt06x/ANwxmydgiF83mTEvMMY9cc4JLNWBJsD/FzmylSHPKAbC8QWutLSqlxwAUgHlgL7AWitdYppmYRQPHsjk0Ikb7JeyaTlJrEO8HvYGdjh18BP+Y8P4eLv0TS/Y+eeCUNRLusIC7lFD+2/JG+ga9bOmSRBZTWOns3qFR+YDHQBYgGFpqmP9ZalzO18QH+1Fr7p7N8X6AvgLe3d+D8+fMzFUdMTAyurq6Pb5jHWGPe1pgzZF3eSYYkuu7qip+rH2P9xxrnpWpWhyez5PQtrjgOI1GF42zrzMeVPia4QPBTbzOzZF9nTKNGjfZqrR8/sJXWOlt/gM7A9Pumu2PsSooC7Ezz6gJrHreuwMBAnVmbNm3K9LK5mTXmbY05a511ec8+MFszGr3m9Bodm5isp245o2t9vk6XGrpC95uzRx+8dEb3WdpH77+yP0u29zRkX2cMsEdn4PvaEucgLgB1lFL5MHYxNQH2AJuAThivZOoBLLVAbEKI+2it+W7nd1QqWInIa+V5du4mbsYmUbdMAf7TpTrPlCsIwLR20ywcqTAHS5yD2KWUWoTxUtYUYD8wBVgJzFdKfW6aNz27YxNCPGjNmTXsv7qfaq7vM2rpUeqU8eL95hUILOVl6dBENrDIVUxa64+Bjx+afRaobYFwhBDpWHxkJV2XdMbOUAzbuPr8p0t12lcvhlJys5u1kDuphRAPSEhOZeDSyUw9Mgh77cPb/jP4qFU9PPLZWzo0kc2kQAghAEhJNbBobwSj183mRPLHFHaqyp/dVhDo42Pp0ISFSIEQwspprVl95CpfrznB2ai73HKdga+nH0f6b8fFwcXS4QkLkgcGCWHFwq7coevUnfT7bR8Otjb0bXabu6nn+Kzxh1IchBxBCJHXzdg/g92XdlOhQAXKFyhPxYIV8XIswXfrTzFn53k8nO35vENVutYuSePZIfi4+9ClShdLhy1yACkQQuRh0QnRDPhzACmGFJINyWnzFY7YG0oQXKoJ/+32HQXd8rEzYidbzm/hu+bfYW8rJ6RFBruYlFLllVIblFJHTNPVlFIfmjc0IcTT+u3Qb8SnxLO9z3bWvHiSZz1+xCvpbco6t6d68WL8FTmFLkvacjP+Jt9s/4b8Tvl5reZrlg5b5BAZPYKYCrwP/AygtT6klJoLfG6uwIQQT0drzZR9U6hSqDo/r4d1x07i7e7H5I5teKFmCWxsFLMOzKLvir4ETgnkfPR5RtYfiauD9Y1lJNKX0QKRT2u9+6EbZFIe1VgIYXlhd8I4FHmIQqlvsUPd4P3mFehdrzTODrZpbXpU70H5AuXp+HtHHO0ceTv4bQtGLHKajBaIKKVUWYxPgEMp1QnjMxyEEDnQ1dsJjDm8FKUdec73BcZ3rkthN6d029b1qcuhfoe4FnuNwi6FszlSkZNltEC8hXG8pIpKqUvAOeAVs0UlhMgUg0Ez7+8LfLFqL5fYTAOfDvzaK+Sxw2MUdiksxUH8Q4YKhNb6LNBUKeUC2Git75o3LCHEkzpx9S4jlhxm7/lbFPLeib6TyNctBsvYSSLTMnoV05dKKU+tdazW+q5SKr9p1FUhhIVdvZ3AsMWHaDlhC6euRdI0aD+XDHMo61KWWsVqWTo8kYtltIuppdZ6xL0JrfUtpVQrQC51FcJCIu/cpeEvbTh/6ypK58PH25OzcXs4fDSGWsVq0dO7pxw9iKeS0aE2bJVSjvcmlFLOgOO/tBdCmElsYgo/bDpNtfHvcCJ6CwVdPKhUzAE7h2u8UPl5dr22i92v76aye2VLhypyuYweQcwBNiilfsF4JVNvYJbZohJC/IPWmj8PX+XjZUe5GnOJm/nm07RUG9b1WG7p0EQeldGT1F8rpQ5jfDyoAj7TWq8xa2RCiDTX7iTw4R9HWHssEv/iHhQqtYKo85op7b63dGgiD8vwWExa61XAKjPGIoR4SKpBM3f3Bb5efZykFAPDW1akYqnLhMxaxIf1P6R0/tKWDlHkYRkqEEqp54GvgMIYjyAUoLXW7maMTQirdjjiNh/+cZj9EdcoU+w8LQIcuGV7kH4rZ1PCvQTDnh1m6RBFHpfRI4ivgbZa6zBzBiOEgPikVMavO8Gk7UtJdgwl3m0HF2/dZXOo8XNPJ09md5gtz2sQZpfRAhEpxUEI89t19gZDFx/iSPRybjh8h7u9O10qd6JLlS5UKlSJgvkKks8+n6XDFFYiowVij1Lqd+APIPHeTK31f80SlRBWJjYxha9WH2f2jvMUyZ9AkusM6nnXY3339TjZpT+GkhDmltEC4Q7EAc/dN08DUiCEeErbTkcxdPEhLkXH0/MZXw7EjeLA6QSmt5suxUFYVEYvc+1l7kCEsDY3Y5MY82cYC/dGUKagCwveqMv5uI18snAJY5qMoULBCpYOUVi5jF7F5AT0AaoAaX/SaK17mykuIfIsg0GzaG8EX64KIyYhhTcblmVQUz/iUm7T5oe3qFm0Ju89856lwxQiw11MvwLHgebAp0A3QE5aC/GETly9y4d/HObv8FvU8s3PFx39Ke/tBsCQtaO4EXeDNa+swc5GHhcvLC+j/wrLaa07K6Xaa61nmR43KndSC5FBcUkpfL/hNNP+Oourkx1fveBP50AfbGyMg+kdijzE5L2T6R/Un+pFqls4WiGMMlogkk2v0UqpqsBVwDezG1VKeQLTgKr8f2ynE8DvpvWGAy9qrW9ldhtC5ASpBs3ivRGMW3uCa3cTeTGoBMNaVsLLxSGtjdaagasHkt8pP580+sSC0QrxoIyO5jpFKZUf4/Dey4BjGO+szqwJwGqtdUUgAGN31TBgg9baD9hgmhYix/th9w+M2jiKa7HXHpi//XQUrb//iw8WH6KYpzOL3qzL150CiIg5zvtr3+d41HEAFoctJjQ8lM8afYaXs5clUhAiXRk9gthg+mt+C1AGQCmVqUFglFLuQAOgJ4DWOglIUkq1B0JMzWYBocDQzGxDiOwy7/A8BqwaAMC3O76lb2Bfmvi2Y/b2s2w5HUlJ91JMerkxrf2LopTiTuIdOszvwLnoc4zbMY72Fdqz/+p+qnlXo29gXwtnI8SDMlogFgM1H5q3CAjMxDbLANeBX5RSAcBeYCDgrbW+AqC1vqKUkgfkihxt7+W99F7Wm/ol6/Nj6x/5ets4Ju6exIRdE4wNHOF2iiM4LUKpNgD0X9mfC7cvsOylZey5vIdJf0/iZvxNZnWYha2NrQWzEeKflNb60R8qVRHjpa1fA+/f95E78L7WusoTb1CpIGAnUE9rvUspNQG4A7yttfa8r90trXX+dJbvC/QF8Pb2Dpw/f/6ThgBATEwMrq6umVo2N7PGvM2R882km7y5700Uih+q/8SxKDeWnk4mMjGSoh4XaV7ShYL5bJl8djKnY04zouIIUnQKY46PoZdvL7qX6g5AfGo8l+MvU9a1bJbGB7KvrcmT5t2oUaO9Wuugx7V7XIFoD3QA2mE893DPXWC+1np7hiP6/zqLADu11r6m6foYzzeUA0JMRw9FgVCt9b/eKRQUFKT37NnzpCEAEBoaSkhISKaWzc2sMe+szjkpNYnGsxqz78o+xjdaxsIddpy5HktACQ/eb16ReuUKpD3q807iHdrOa8tf5//Cyc6JWsVrsbH7xmw5WpB9bT2eNG+lVIYKxL92MWmtlwJLlVJ1tdY7Mrz1f1/nVaXURaVUBa31CYwPITpm+ukBjDW9Ls2K7QmR1d5d/S7bLm7jmfyfMnZZImUK2jH5lUCaV/H+xzOg3R3dWdVtFS8ufJHdl3Yzp+Mc6UoSuUZGz0F0VEodBeKB1RivPBqktZ6Tye2+DfymlHIAzgK9MF5RtUAp1Qe4AHTO5LqFMJuf90zjxz0/4pHyPHHRwYxu60e3OqWwt330BYH57POxvOtyElMTZWwlkatktEA8p7X+QCnVEYjA+OW9CeOzqp+Y1voAkN7hTZPMrE8Ic0hISaDr4q7Y2djRsFRD4uJdGR7aH6fUAPr4j2BUG3888zk8fkWAUkqKg8h1Mlog7E2vrYB5WuubDx9KC5HXjNwwkj+O/0FR1+IsOrYIAEe8+b3zfNr5V7RwdEKYX0YLxHKl1HGMXUz9lVKFgATzhSWEZW08t5HxO8dTr8jLXI/ohqOOpE7FKIY3bUu1IlIchHXI6HDfw5RSXwF3tNapSqlYoL15QxPCfFIMKSw4uoDWfq3xcPJ44LNb8bfosvBVnCjBhXPP09a/CMNahuDjJU9yE9blXwuEUqqx1nqjUur5++bd30QeGCRypTmH5tBraS+qFq7Kny//iY+HDwaDZs2xS/Rd0YuohKvUdfuR714NIbhMAUuHK4RFPO4IogGwEWiLcVA99dCrFAiR62it+X7X9/i4+3Dh9gXqTK/Dm1WmMG/fDk4lzCDF5gqd/AYzv2tfbG3kXJuwXo8rEHeVUoOBI/y/MGB6L0SutP3idvZf3c9PrX8iv50/PZd15KOdxqEwfD0r8Z+Wk2lXoe0/7mkQwto8rkDcu3e7AlAL481rCuMRxRYzxiWE2UzcPRF3Rw/Onq/Bwj1RVHabiFexxfQK7MBLVV/CRmV0kGMh8rbH3Un9CYBSai1QU2t91zQ9Glho9uiEyGKnb5xn4bHF5De0Y9GeKHo9U5rBz5XH1bGLpUMTIsfJ6GWuJYGk+6aTeIoHBgmR3bTWLDt4mX7LRmIwpNLU51U+a9MAP9PjPoUQ//Qkz6TerZRagvH8Q0eMz2wQIltdjbnKrfhbVCpU6ZFtbifcZuuFrbQu3xqA63EGevzyN5tPXuZGvpXUL9mc+b07ZFfIQuRaGeps1Vp/gXG8pFtANNBLaz3GnIEJ8bC45Dga/NKAGj/XYMv5R58CG7d9HG3mtWHnxd1M2XKGkdvi2Rt+k+ZBZ0nS0YxqODgboxYi98rw2Tit9T6t9QTTz35zBiVEeoauG8qpm6co7FKYtvPacvDqwXTb/Xn6TwA6z/mUL/88TiUvW9YNbsih6AWUL1CeJmVkyC8hMkIu1xC5wvqz65n09yQGBg9ka++tuDu603xOc87cPPNAu+PXLrDvyj6UduZy0nq+6VyWQTUduRJ3jF2XdtE/qL9cpSREBsn/FJHjRSdE02tpLyoUqMCYJmMo6VGSta+sJdmQTJt5bUg1pJJq0Py68zzP/TQRgDYlh2EgkSspa1FK8cPfP5DPPh89qvewcDZC5B5SIESOlmJI4bVlr3Hl7hVmd5yNs70zAJUKVeLHVj9yPOo4E7cvpM3ErYz64wjK6SAFnQvzR68PqV28NpP3TOZ28m3mHZnHK/6v4Onk+ZgtCiHukQIhcqwUQwo9/ujB4rDFfN3sa2oXr/3A5w18WpPPtiDD144jOi6J77sGEKP20rp8S2yUDf2C+hEWFcbXJ74mISWBt2q/ZaFMhMidpECIHOlecZh7eC5jmoxhcN0HrzxafyySNt/vwDGpGYm2+5jauwSFvC5yM/4mLcu1BKBLlS54Onmy/cZ2ni35LNW8q1kiFSFyLSkQIkd6d/W7acVh2LPD0ubHJaUwYslhXpvu6yE9AAAfFklEQVS9h8LuTizt8TE2yoZZB6ex6vQqbJQNzco2A8DZ3pmeAT0BeKuWHD0I8aQyeqOcENnGoA3MOTyHbv7dHigOW09FMWLJYS7eiuONBmUY/Fx5HO1s6VCxAzMOzKCEewmCiwfj5eyVtszw+sOJuRbDC5VesEQqQuRqcgQhcpyTN04SnRBN49KNAYiOS+K9hQd5ZfoubG0U816vw/BWlXC0swWMRwc3429yKPJQWvfSPYVdCtOtZDfsbe3/sR0hxL+TIwiR4+yK2AVAcPFgFu65yJhVx7kTn8xbjcrydmM/nOxtH2gf4htCpYKVCIsKo6Vfy/RWKYTIBCkQwqIm75lM3RJ1CSgSkDZv16VduDq4M3pxNH+fDyewVH4+71CVSkXd012HUopPG33KjP0zqFm0ZnaFLkSeJ11Mwuxuxd8iZGbIP8ZPOnD1AP1W9uPj0I//3zY2iSVHQ0lNKM2p67F89YI/C9+o+8jicE+nyp34s9ufcpe0EFlI/jcJsxu7dSybz2/m3TXvovX/H0b47Y5vAVh7Zi2xiXH8uiOcBuNWczXuJDWL1GbjkBC61CqJjTz2UwiLkAIhzOrSnUt8v/t7SnmUYt+VfSw7sQyAi7cvMv/IfGoWrUl8Sjwtf/6RUUuPUqTAFVCpfNC4DV4uDhaOXgjrJgVCmNUnmz8h1ZDK+u7r8fPy4+PQjzFoAxN2TUBrzUtlv8VGO3Pk5ka+fqEazWveBYwnqIUQliUFQpjNiagTzNg/g35B/SjnVY6PG37MwciDzDwwk8l7plDEvhE/rIvFJ98z2Lvso1NQcXZf2k0pj1J4u3pbOnwhrJ4UCGE2IzeOxNnemZENRgLwUtWX8POqwGvL+hKbfBeP1I5MeKk6n7foybW4SPZc3sOuS7sILiFHD0LkBBYrEEopW6XUfqXUCtN0aaXULqXUKaXU70op6YDOxY5cO8LisMW8V/c9CrsUJtWg+W3XRRJuPI8mlTJuwex6vy/tqxenlV8rbJUtP+/5mQu3L1CneB1Lhy+EwLL3QQwEwoB71y9+BXyntZ6vlJoM9AF+slRw4ulM3zcdB1sHBtQewOGI2wz77yGOXr5DvTJtcC0cR7/g7rg6Gv/5eTl7Ub9UfWYdND7mXI4ghMgZLFIglFIlgNbAF8BgpZQCGgMvm5rMAkYjBSJXSkpNYs7hObT2a8uU0GtM/essBV0dmfRyDVr7F0WpZ/6xTLvy7QgND8XOxo4aRWpYIGohxMMs1cX0H+ADwGCaLgBEa61TTNMRQHFLBCae3oqTK4iKi+LUuSB+3nKWF4N8WDe4IW2qFcP4t8A/tavQDoAA74C0hwIJISwr248glFJtgGta671KqZB7s9NpqtOZh1KqL9AXwNvbm9DQ0EzFERMTk+llc7OsyHv3zd38Ev4LQysMxdfF94HPrscZeGffV9hqL5wTq9CvlhOVC9xk/65tj11vsFcwAU4BWb5fZF9bD2vMGcyYt9Y6W3+AMRiPEMKBq0Ac8BsQBdiZ2tQF1jxuXYGBgTqzNm3alOllc7OnzTsuKU6X+q6UZjS6yLgi+kTUCa211rGJyfrbtSd02Q/naD620c2mv6ETklOyIOKnJ/vaelhjzlo/ed7AHp2B7+ts72LSWg/XWpfQWvsCLwEbtdbdgE1AJ1OzHsDS7I5NPN647eM4f/s8P7b6kVRDKo1nNeaHLVtpNC6U7zecorD3blAGJrUfnDYctxAid8pJ90EMxXjC+jTGcxLTLRyPeMjF2xcZs3UML1R6gb6BbzIyeC7XYu4ycGN7btj9wsA2t7ml1lLPpx7lC5S3dLhCiKdk0eG+tdahQKjp/Vmg9r+1F5Y1dP1QNJqXyo+g2XebOXs9kcCC47njNJVTtxYyeMNcAKa1nWbhSIUQWUGeByEyZPvF7cw7Mo9Gxfrx/u9XKF3AhR+71aRFlSLY2PQhPjme3Zd2c/LGSV4NeNXS4QohsoAUCJEh3++cip1y4fSZJnSrXZKP2lTG2eH/5xic7Z1p6NuQhr4NLRilECIrSYEQ/yopxcDPm0+z+Ogq8hHAz92eoaV/UUuHJYTIBlIgxCPtCb/JyCVHOHLtJClOkQxvPEyKgxBWJCddxSRyiEvR8bw9bz+dJu/gTkIynZ65CcDzlVtZODIhRHaSIwiRJiE5lZ9CzzB58xkA3mnix5sNy/DqHz9S0qMkfl5+Fo5QCJGdpEAItNasD7vGJ8uPEnErnrYBxRjWsiLFPZ1JNaSy8dxGXqj0wiPHURJC5E1SIPK4dWfW4eHkQe3i6d9icvFmHKOXHWXD8Wv4FXZl3ut1qFu2QNrne6/sJTohmqZlmmZXyEKIHELOQeRh0QnRtJ7bmuBpwXRZ1IWzt86mfZacamDy5jM8990Wdpy9wchWlfjvW7X4/dRndFnUhcSURMBYYACalG5ikRyEEJYjRxB52PITy0k2JNOrei9+P/o7S8KW0MK7BVH57Zm2KZXjV+/StJI3n7Svwt2U8zz7S10ORR4CwMXehentprP+3HpqFKlBIZdCFs5GCJHdpEDkYYvCFuHj7sP0dtP5vPHnjNr4MbMOzGb5H8txU/60rtGY/Pnt+XxrDL8e+hVne2f+fPlPdkTs4LMtn1HaszTbLmxjUJ1Blk5FCGEBUiDyqDuJd1hzeg39gvqhlOLCdUdOnXiZovGNKV9mD+fil7LgxESc7JxwsnOifqn6TGs7jeLuxWlerjnHrh/jo9CPAGhWppmFsxFCWIIUiDxq5cmVJKYm0qZ8B0YvO8rM7eGUKpCPkbUL8+bz3wLforVO98okG2XDrA6zOHPrDCdvnOTZks9mfwJCCIuTApFHLTq2CC8nbz5akMyl6HB6PuPLBy0qsHv71rQ2/3bZqouDCxu7b+T87fPyCFAhrJQUiFzs2PVjVChQAVubBx/Ms+NsBH8cX4FLSjPcCjiw4I2a1C7t9cTrz++cn/zO+bMqXCFELiOXueZScw7NocqPVei/sn/avDsJyXz4x2HaT/8BA0kMeuZVVr5TP1PFQQgh5AgiFzoRdYI3V7xJfqf8TNk3hZpFa+Lj2I7Ry44SFZNI4aIH0ImF+Lh5J2xt5O5nIUTmyBFELhOfHM+Li17Eyc6JA28eoE6xRvRbOYA+c+fi4ZJM4+ANnLy9gRcqPf+PrichhHgScgSRywxZO4RDkYeY0nIhXyy7RsSZvjg4HyfOfSz7kzTrD97gZf+X+azRZ5YOVQiRy0mByEU2ndvET3t+IsirB2OWOJPP4QbvNwskuMIKms9pTPUiQXzV9CsCiwVaOlQhRB4gBSKHSDGkEHEnAl9P33Q/T0hO5bUlQ7HVXkRHduS1+mV4s2FZvFwcAIj6IAo7G9mdQoisI+cgcohvtn1D+YnlCY8O/8dnoSeuETxuImfv/s2z3n346/3mjGhVKa04AFIchBBZTr5VcgCtNdP2TyPZkMyUvVP4ssmXAETciuPzFWGsPnqVO65zyO9UiD9fG00+eyfLBiwEkJycTEREBAkJCZYOJY2HhwdhYWGWDiPbPSpvJycnSpQogb29fabWKwUiB9h6YStnb50lv1N+pu+fzrB6o5i5LYIfNp1GKXihTizjD+7lm8bfkM8+n6XDFQKAiIgI3Nzc8PX1zTEPk7p79y5ubm6WDiPbpZe31pobN24QERFB6dKlM7Ve6WLKAX458AuuDq5MbTuVa7HXqPufMYxfd5KmlbzZMCSEsNiZFHAuwJtBb1o6VCHSJCQkUKBAgRxTHMSDlFIUKFDgqY7w5AjCwmKSYlhwdAEdKnRm6+HS2Bm8uZS8jOV93qa+XyFWn17NqtOr+LLxl7g6uFo6XCEeIMUhZ3va/SNHENnswNUDTNk7BYM2ALD42GJik2PZfbQayw9epVmpbtw2HKSgZxSbwzfz/O/P41/YnwG1B1g4ciFyt5kzZzJgQOb+H/Xs2ZNFixY9dv2XL19Om37ttdc4duxYpraXU2R7gVBK+SilNimlwpRSR5VSA03zvZRS65RSp0yvuXqUuFRDKsmpyQ/M01rTfUl33ljxBh1/78j+i5cYunoidoaiVPKqxYp3nmVWl2HY29gzaM0gWs9tja+nL+u7r8fN0fr6VYXITR4uENOmTaNy5coWjOjpWeIIIgUYorWuBNQB3lJKVQaGARu01n7ABtN0rrTuzDrKfl+WkFkhpBpS0+avPLWSw9cO82zxViw/sZLa04KJTNxL23Jd+W//elQs4k4hl0J0qtyJtWfWUsytGBu6b6CwS2ELZiNEzjZ79myqVatGQEAAr7/+OsuXLyc4OJgaNWrQtGlTIiMj/7FMZGQkHTt2JCAggICAALZv3054eDhVq1ZNazNu3DhGjx79j2U//fRTatWqRdWqVenbty9aaxYtWsSePXvo1q0b1atXJz4+npCQEPbs2QPAvHnz8Pf3p2rVqgwdOjRtXa6urowcOZKAgADq1KmTbqyWlO3nILTWV4Arpvd3lVJhQHGgPRBiajYLCAWGprOKHOtW/C2GrB3CLwd+oYhrEbZf3M7Pe3+mf63+aK0Zuf4znG2KcOF0X8q7NuSi3RhSUxT/aTfwgUH1RjUYhUEb+KbZNxR1K2rBjITImE+WH+XY5TtZus7Kxdz5uG2Vf21z9OhRvvjiC7Zt20bBggU5f/487u7u7Ny5E6UU06ZN4+uvv+bbb799YLl33nmHhg0bsmTJElJTU4mJieHWrVsZimvAgAF89JHxaYuvvvoqK1asoFOnTkyaNIlx48YRFBT0QPvLly8zdOhQ9u7dS/78+Xnuuef4448/6NChA7GxsdSpU4cvvviCDz74gKlTp/Lhhx8+wW/JvCx6DkIp5QvUAHYB3qbica+I5Ko/mzee24j/T/7MPjib4c8O5+w7Z2lSugkjNozgyJXzvDhzBoeu76YwnfjuxUAOD3+PsAEHWd99PSU9Sj6wrkqFKjG/03x8PHwslI0QucPGjRvp1KkTBQsWBMDLy4uIiAiaN2+Ov78/33zzDUePHk13uX79+gFga2uLh4dHhre5adMmgoOD8ff3Z+PGjemu/35///03ISEhFCpUCDs7O7p168aWLVsAcHBwoE2bNgAEBgYSHh6e4Tiyg8WuYlJKuQKLgUFa6zsZPduulOoL9AXw9vYmNDQ0U9uPiYnJ9LL3SzIkMePcDH6P+B0fZx9+qPEDFWwrsGvbLl4t8Cqh57ZQ76deJHMHZ3sPJtRuhced02z76zQANtgQev7p48iorMo7N7HGnMH8eXt4eHD37l0ABoeUfEzrzLm3/keJj48nKSkprV1qair9+/dnwIABtGrVir/++osxY8Zw9+5dEhIS0tpqrbl79y5JSUlp60pISCAlJSVtXbdv306bTk5OJj4+nuvXr9OvXz82b95MiRIl+PLLL7l9+zZ3794lNTWV2NjYB2KJjY0lLi6O5OTktPn3x2Fvb09MTAwASUlJxMfHPzbn9KSmpj5yuYSEhEz/O7BIgVBK2WMsDr9prf9rmh2plCqqtb6ilCoKXEtvWa31FGAKQFBQkA4JCclUDKGhoWR22fu1+q0VqyJW0S+oH+OeG5d2I9vBi9GMX2qHa/Lz3Lb/HYAvQr6gff2WT73Np5FVeecm1pgzmD/vsLAwi9+U1rp1azp27MiwYcMoUKAA58+fJyYmhnLlyuHm5sbChQuxtbXFzc0NJycnHBwccHNzo2nTpsyZM4dBgwalfZGXLVuWqKgokpKScHV1Zd26dbRo0QI3Nzfs7e1xdnbG3t4epRS+vr6kpqayfPlyOnXqhJubG56enhgMhrTfia2tLS4uLoSEhDBs2DASExPJnz8/S5Ys4e23305rd+/13voz8zv9txsEnZycqFGjRqZ+v5a4ikkB04EwrfX4+z5aBvQwve8BLM3u2J5UfHI8a86sYXCdwfzY+kfy2efjUEQ0fWb+TfsftnE5OoGpz39Bac/SuDm40b9W/8evVAiRYVWqVGHkyJE0bNiQgIAARowYwejRo+ncuTP169dP63p62IQJE9i0aRP+/v4EBgZy9OhR7O3t+eijjwgODqZNmzZUrFjxH8t5enry+uuv4+/vT4cOHahVq1baZz179uTNN99MO0l9T9GiRRkzZgyNGjUiICCAmjVr0r59+6z/ZZiD1jpbf4BnAQ0cAg6YfloBBTBevXTK9Or1uHUFBgbqzNq0aVOml71nz6U9mtHoRUcX6VORd3WfmX/rUkNX6IBP1uiJG07qO/FJWmutw2+F632X9z319rJCVuSd21hjzlqbP+9jx46Zdf2ZcefOHUuHYBH/lnd6+wnYozPwfW2Jq5i2Ao864dAkO2N5WgcjDwKw8bATHxzcQj57W957rjw9nvHFzen/g2OV8ixFKc9SlgpTCCEyRYbayCSDQbPgwFaUdmTVQQOvBJdmYBM/Crg6Wjo0IYTIElIgMuHY5TuMWHKY0Mi9eDmVZW3/EMoVljudhRB5ixSIJ3A3IZnx604ye8d5PJzscHC+SMeqL0hxEELkSVIgMuBOwh3KTKiIW9IrENeQl2uX5JV6blT+6RYB3gGWDk8IIcxCCsRjhF25Q8/fx3Ej4QrKbgPr+48gwMeT1adXA1DNu5qFIxRCCPOQ4b4fITYxhdHLjtJm4lbCbv8JQLThEL6FjBdgHYo8BIB/YX+LxSiEEOHh4SxYsMAs65YCkY7d527ScsJfzNoRTtsa+YhTBwjxDSHFkJJ25HAo8hA+7j7kd87Vo5ILIbJQSkpKtm8zPDychQsXmmXdUiDuk5Ccyhcrj9Flyg4AFrxRl7IlD5OqU5nQYgIF8xVk+cnlgLFASPeSEJYVHh5OxYoVee2116hatSp9+vRh/fr11KtXDz8/P3bv3k1sbCy9e/emVq1a1KhRg6VLl6YtW79+fWrWrEnNmjXZvn07AFeuXKFBgwZUr16dqlWr8tdffwHGobnvWbRoET179gSMd1APHjyYRo0aMXTo0Edub+bMmXTo0IG2bdtSunRpJk2axPjx46lRowZ16tTh5s2bAJw5c4YWLVoQGBhI/fr1OX78eNp23nnnHZ555hnKlCmT9gCjYcOGsWPHDqpXr853332Xpb9fqz4HYdAGvt3+LZ//9TkKO5KS8mFI8aRb1Q/5vlNzXBzteGfDXKp5V6OadzVa+bVi+YnlxCXHERYVRpvybSydghA5wqDVgzhw9UCWrrN6ker8p8V/Htvu9OnTLFy4kClTphAYGMjcuXPZunUry5Yt48svv6Ry5co0btyYGTNmEB0dTe3atWnatCmFCxdm3bp1ODk5cerUKbp27cqePXuYO3cuzZs3Z+TIkaSmphIXF/fYGE6ePMn69euxtbVlxIgR6W4P4MiRI+zfv5+EhATKlSvHV199xf79+3n33XeZPXs2gwYNom/fvkyePBk/Pz927dpF//792bhxI2AsXlu3buX48eO0a9eOTp06MXbsWMaOHcvq1auf7heeDqstENcTr9Ps12ZsPLeRsm71ibzljqtjHAbXY6y48g4jYoMgFnZG7OSrpl8B0LZ8W2YfnM2M/TNIMaTIEYQQOUDp0qXx9zeeC6xYsSJNmjRBKYW/vz/h4eFERESwbNkyxo0bBxhHN71w4QLFihVjwIABHDhwAFtbW06ePAlArVq16N27N8nJyXTo0IHq1as/NobOnTtja2sLwNq1a9PdHkCjRo1wc3PDzc0NDw8P2rZtC4C/vz+HDh0iJiaG7du307lz57R1JyYmpr3v0KEDNjY2VK5cOVseLmSVBWLdmXX0/LsPSYYUCiYPJPV6U4Y9W4bBzSpw4c5p6s2oR/M5zWnl1wqArlW7AvBc2eewt7Hn621fA3IFkxD3ZOQvfXNxdPz/6AU2NjZp0zY2NqSkpGBra8vixYupUKHCA8uNHj0ab29vDh48iMFgwMnJCYAGDRqwZcsWVq5cyauvvsr7779P9+7duf+RBAkJCQ+sy8XFJe291jrd7e3ateuxsRoMBjw9PTlwIP2jsfuXNw6pZF5WeQ7iyAU7UpN8Kat/4N1n3mDL+40Z2boyzg62VChYgRUvr+Dy3ctM3D2RBqUapD24x93RnRDfEC7euYiDrQPlC5S3cCZCiMdp3rw5EydOTPtC3b9/P2B83kPRokWxsbHh119/JTXV+Hjg8+fPU7hwYV5//XX69OnDvn37AOPzZ8LCwjAYDCxZsuSJt5cR7u7ulC5dOu2ks9aagwcP/usybm5uac+UyGpWWSB6Bz/LML9v2T+8OyNaVcLHK98Dn9cpUYeFnRfiYOvAG4FvPPBZ2/LGQ8IqhapgZ2OVB2BC5CqjRo0iOTmZatWqUbVqVUaNGgVA//79mTVrFnXq1OHkyZNpRwGhoaFUr16dGjVqsHjxYgYOHAjA2LFjadOmDY0bN6Zo0Uc/CvhR28uo3377jenTpxMQEECVKlXSTnI/SrVq1bCzsyMgICDLT1Jn+3DfWflj7uG+45Li/jHv3K1zmtHoHkt6ZHrblmSNQ19bY85ay3Df1iTPDPedmzjbO/9jnq+nL181/YpGvo0sEJEQQmQfKRCZ8EG9DywdghBCmJ1VnoMQQgjxeFIghBCZprPhUkuReU+7f6RACCEyxcnJiRs3bkiRyKG01ty4cSPt/o7MkHMQQohMKVGiBBEREVy/ft3SoaRJSEh4qi/E3OpReTs5OVGiRIlMr1cKhBAiU+zt7SldurSlw3hAaGgoNWrUsHQY2c5ceUsXkxBCiHRJgRBCCJEuKRBCCCHSpXLzFQhKqevA+UwuXhCIysJwcgtrzNsacwbrzNsac4Ynz7uU1rrQ4xrl6gLxNJRSe7TWQZaOI7tZY97WmDNYZ97WmDOYL2/pYhJCCJEuKRBCCCHSZc0FYoqlA7AQa8zbGnMG68zbGnMGM+VttecghBBC/DtrPoIQQgjxL6yyQCilWiilTiilTiulhlk6HnNQSvkopTYppcKUUkeVUgNN872UUuuUUqdMr/ktHas5KKVslVL7lVIrTNOllVK7THn/rpRysHSMWUkp5amUWqSUOm7a53WtYV8rpd41/fs+opSap5Ryyov7Wik1Qyl1TSl15L556e5fZfS96fvtkFKqZma3a3UFQillC/wAtAQqA12VUpUtG5VZpABDtNaVgDrAW6Y8hwEbtNZ+wAbTdF40EAi7b/or4DtT3reAPhaJynwmAKu11hWBAIy55+l9rZQqDrwDBGmtqwK2wEvkzX09E2jx0LxH7d+WgJ/ppy/wU2Y3anUFAqgNnNZan9VaJwHzgfYWjinLaa2vaK33md7fxfiFURxjrrNMzWYBHSwTofkopUoArYFppmkFNAYWmZrkqbyVUu5AA2A6gNY6SWsdjRXsa4wDjjorpeyAfMAV8uC+1lpvAW4+NPtR+7c9MNv0+OmdgKdSqmhmtmuNBaI4cPG+6QjTvDxLKeUL1AB2Ad5a6ytgLCJAYctFZjb/AT4ADKbpAkC01jrFNJ3X9nkZ4Drwi6lbbZpSyoU8vq+11peAccAFjIXhNrCXvL2v7/eo/Ztl33HWWCBUOvPy7KVcSilXYDEwSGt9x9LxmJtSqg1wTWu99/7Z6TTNS/vcDqgJ/KS1rgHEkse6k9Jj6nNvD5QGigEuGLtXHpaX9nVGZNm/d2ssEBGAz33TJYDLForFrJRS9hiLw29a6/+aZkfeO9w0vV6zVHxmUg9op5QKx9h92BjjEYWnqRsC8t4+jwAitNa7TNOLMBaMvL6vmwLntNbXtdbJwH+BZ8jb+/p+j9q/WfYdZ40F4m/Az3SlgwPGk1rLLBxTljP1u08HwrTW4+/7aBnQw/S+B7A0u2MzJ631cK11Ca21L8Z9u1Fr3Q3YBHQyNctTeWutrwIXlVIVTLOaAMfI4/saY9dSHaVUPtO/93t559l9/ZBH7d9lQHfT1Ux1gNv3uqKelFXeKKeUaoXxr0pbYIbW+gsLh5TllFLPAn8Bh/l/X/wIjOchFgAlMf4H66y1fvjkV56glAoB3tNat1FKlcF4ROEF7Ade0VonWjK+rKSUqo7xpLwDcBbohfEPwDy9r5VSnwBdMF61tx94DWN/e57a10qpeUAIxlFbI4GPgT9IZ/+aiuUkjFc9xQG9tNZ7MrVdaywQQgghHs8au5iEEEJkgBQIIYQQ6ZICIYQQIl1SIIQQQqRLCoQQQoh0SYEQ4gmYRk3tb3pfTCm16HHLCJFbyWWuQjwB07hWK0yjhwqRp9k9vokQ4j5jgbJKqQPAKaCS1rqqUqonxtE0bYGqwLcYb1p7FUgEWpluYvpfe/erEkEUxXH8eywmo09gFV2w+R4WBYsYbQq+gMlu8BWMPoBJNosIFvEJDLJg3GO4F1nkwu4Mw4Tl+2nzlzth+HEvM+fsUMrNb1N+YjrPzPfxH0NaziUmqZtr4CMzJ8DVv2O7wDGlpPwN8FOL502B03rOPXCRmQfAJXA3yqilHpxBSMN5qr03ZhHxDTzW/a/AXq2sewg8lGoIAGyOP0xpNQaENJzFej/zhe055V3boPQqmIw9MKkPl5ikbmbAVp8Laz+Oz4g4gr/ewftDDk4akgEhdZCZX8BzbR5/2+MWJ8BZRLwAb6xhu1utDz9zlSQ1OYOQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqekXKVgH4Pp2BjEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10879ac88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd8jef/x/HXlSVBltEYUZsiEsTqImZRo7VLa/ysfu2qFvVtqVapUbRKBy06bLVaSomqUZq0dmxBjBgxMkTW9fvjHL7KISeRkzs55/N8PM4jOfe57/t8rtyc97nucd1Ka40QQghxPyejCxBCCJEzSUAIIYSwSAJCCCGERRIQQgghLJKAEEIIYZEEhBBCCIskIIQQQlgkASGEEMIiCQghhBAWuRhdwOMoVKiQLlWqVKaWjY+PJ1++fFlbUC7giO12xDaDY7bbEdsMGW93eHj4Fa114fTmy9UBUapUKcLCwjK17JYtWwgJCcnagnIBR2y3I7YZHLPdjthmyHi7lVKnrZlPdjEJIYSwSAJCCCGERRIQQgghLJKAEEIIYZEEhBBCCIskIIQQQlgkASGEEMIiCQghhMhNtIa0tGx5KwkIIYTITdauhcBAiIy0+VtJQAghRG6hNYwfj05IAH9/m7+dBIQQQuQWmzfDrl3MrNmWHaev2/ztJCCEECIXuJ2Sypk3RxOdvwCLKjckJVXb/D0lIIQQIocLP32N/7y7CP+9u/nz5R78/HZj6lVIdzDWx5arR3MVQgh7duNWMpN/PcwPu85QxKsIu9f+Tpv61SGvW7a8vwSEEELkQL8dimbUT/vJc+4s38TuodbXU8jv7pqtNUhACCFEDhKbmMwHaw+xaucJPvp7MS/vWImTUjD6P/DUU9laiwSEEELkEDuOX+GtZfu4fvkaGzZN5sl9u1E9e8LYsVCiRLbXIwEhhBAGi01MZsK6w/y46wylC+Vj+/5v8Nn/FyxYAK++alhdEhBCCGGgnSeuMnzpXqKvxTGoeiH6t62DRwNPONkDOnQwtDYJCCGEMEBqmmbGpmOsW7KJD/9awvOn/sGlTWvo/AwEB5seBpOAEEKIbHbhxi2GLNpDgfVrWLtuOm75PFDt20HbtkaX9i8SEEIIkU201iwLj2Lc2kM0OrCV6SsnQN26sHw5FCtmdHkPkIAQQohsEH0zkXdW7GfT4UvULlWAYb2GQvW8MHw45MljdHkWSUAIIYQNaa1ZGhbFBz8fovS5E2w9tBz/UT/h5O0Fo0cbXd4jSUAIIYSNnI1JYOSKfWw/fpX+Nw4w/MdxOPn4QNRZ8K5idHnpkoAQQogslpqmmb8jksm/HsHZSTHf+yz1pv4XVaMGrFoFRYoYXaJVJCCEECILHb8Ux4jl+wg/fY2QioWZ6n6Ggq8NhDp1YN068PIyukSrSUAIIUQWSE5N46utJ5nx2zHy5nHmk45BvFy9OCqqKLzyCsyaBZ6eRpeZITa7H4RS6hul1CWl1IF7phVQSm1USh0z//Q1T1dKqU+VUseVUvuUUjVsVZcQQmS1A+du0Gbmdib/eoQmlf3Y+EZ92hZMRWltGkPpu+9yXTiAbW8YNA9odt+0kcAmrXV5YJP5OUBzoLz50ReYbcO6hBAiS9xKSmXCLxG0+Xw7l+Nu88WrwXzetQaFb1w2Xd/wxhtGl/hYbLaLSWu9VSlV6r7JbYAQ8+/zgS3ACPP0BVprDfyplPJRShXVWl+wVX1CCPE4dhy/wqif9nP6agKda5VgVPNKeOd1hfh4aN3a9LNvX6PLfCzK9Jlso5WbAmKt1jrA/Py61trnntevaa19lVJrgYla623m6ZuAEVrrMAvr7Iupl4Gfn1/wokWLMlVbXFwc+fPnz9SyuZkjttsR2wyO2e7saHN8smbR4ST+OJeCX15Fjyp5qFTQGQDPQ4co/9lneB49yv6PPiKmTh2b1nJHRtvdoEGDcK11zfTmyykHqZWFaRaTS2v9FfAVQM2aNXVISEim3nDLli1kdtnczBHb7YhtBsdst63b/OvBi4xdeYCY+FRer1+WoY3L4+5qCgeSkqBrV0hNhR9/JLBTJ5vVcT9btTu7AyL6zq4jpVRR4JJ5ehRw790w/IHz2VybEEJYlJicyvifI/juz9NULurFtz1qEVDc2/Tib79BSAi4ucHq1VChQq48IG2JLQ9SW7Ia6G7+vTuw6p7p3cxnM9UFbsjxByFETnDqSjxtZ+3guz9P0+f50qwc8KwpHG7fhsGDoUkT+PJL08zBwXYTDmDDHoRSaiGmA9KFlFJRwBhgIrBEKdULOAPcuRvGL0AL4DiQAPS0VV1CCGGNtDTNgp2RTPr1CG4uTsztXpNGlfwgMREW/AhTp8KhQzB0KPTpY3S5NmHLs5heechLjSzMq4EBtqpFCCEy4sTlOEYs20fY6WvUr1CYCW2rUszHw/Riz56waBEEBpqGzWjd2thibSinHKQWQgjDxSYmMzP0ON9ui8TDzZmpHYJoG1QENW8eNG8OxYubhufu29d03EFZOr/GfkhACCEcXmqaZmnYWaZsOMKVuCTa1fBnRPOKPJFwwzSG0t9/w/jx8M47OeJWoNlFAkII4dD+OXONMasPsi/qBjVL+vJNj1oE+vvApUvQqBGcOgWLF0OHDumvzM5IQAghHNL1hCQ++iWCJWFRFMnvyvwATb06hVHFveHKFWjcGE6ehJ9/hgYNjC7XEBIQQgiHsykimpEr9hN3M56ZSQdpvvoHnCMOgbs7xMWZLnbLkwfWrHHYcAAJCCGEA7lxK5lxaw6x/O8o2sefZMLyj3G9cA6qVoW5c8HXF5ydwc8Pdu0Cp+y+VCxnkYAQQjiETRHRjF6+l7TLlxnYoiaDnwrE9eBPMPdraNbswTOSHDwcQAJCCGHnbiQk8/6ag6wKP8MXobN42jWe/DP+ML24ebOxxeVwEhBCCLsVfjqGwQv3cC0mlvU7Pqd82EZ4/33Q2u6vYcgKEhBCCLuTmqb54vcTfLLxKE1unGDSgZV4bd8M06fDkCFGl5dryE42IYTdcEpKImraLLbVbsrk9YdpHlCEz87+hte2LfDNNxIOGSQ9CCFE7nboECxaRHLEYaqt34BX3A1uP1GS2U38adYoCPX0dPjqSyhSxOhKcx0JCCFE7nLsGMybBx07QlAQ+tgxGD+eiz5F2Fe8Khdf6Ub74d0pm8/NNH/p0oaWm5tJQAghcofUVNMxhNGjISUFSpXiXMkKjLvyBJuHLads8QJ0KJVMr5ceGDBaZJIcgxBC5HwHDkD9+qaRVJs2JfHEKWaWDaHR1C1sjbzBW62qsmbQc5T1cTa6UrsiPQghRM6jNYwbZ7pjm68vLFsGBw6g589nc82mjFsSwemrCTQPKMLoFyvh75vX6IrtkvQghBA5z7hxMHasqecA0Lcvx7fvoVtyRXotCMfFSfFdr9rMfjVYwsGGpAchhMhZvvzSFA49e8Jzz3HjVjJT/4rh+z9P4+nuyphWlXm1bklcneX7ra1JQAghco6ffoL+/eHFF9Fffsnqvef5YG0EMfG3ebVuSd5oXAHfO2cnCZuTgBBC5Ay3bkH//ujatdn6wUw+nxvG7sgYAv29+bZHLar6extdocORgBBCGCsmBry9ScvjztZp85h+SrNn8UGKeLnzQZsqdKlTEmcnGTfJCBIQQojsd+ECrFtnevz6K9GvD6Z/6RaEn06hop8nU5uWoVVQMdxc5DiDkSQghBDZa9Uq6NwZEhNJK16cf+o04d0rfkR7xjO5fSDtavjjJD2GHEECQgiRfbQ2DZMRGEjEuKn8Z28Sp6/dovvTpXijSQW8PVyNrlDcQwJCCGF7Wpvu9ezpSdL87/j692N88vsFini5s7BPXeqWKWh0hcICCQghhG2dOAF9+sC1a+xesJJ3NkZy/FIcbaoVY1ybAOk15GByBEgIYRtpaabB9apWRYeHs6R2Kzp+v5/E5FS+7VGLGZ2rSzjkcNKDEEJkvfh46NYNVqzgUr3G9Kjdk6OuPvSvV4ZBDcvj4SaD6uUGEhBCiKynFMmnIln5ylDeKtGIqv4+rG4XSOViXkZXJjJAAkIIkTW0hmXLSAtpwOJTCUx6cRzx2okRjSvQ5/nSuMjYSbmOIQGhlHoD6A1oYD/QEygKLAIKAH8Dr2mtk4yoTwiRQYcPw6BB8NtvLGvVm1GVX6JumQJMaBtI6UL5jK5OZFK2R7pSqjgwGKiptQ4AnIHOwMfANK11eeAa0Cu7axNCZNDu3dCtGzowkKQ/d/NBs/6Mr/YyH7erysI+dSUccjmjdjG5AB5KqWQgL3ABaAh0Mb8+HxgLzDakOiGEdWbNIm3lSjbXe4kRldpQKagc69oHUszHw+jKRBZQWuvsf1OlhgDjgVvABmAI8KfWupz59RLAOnMP4/5l+wJ9Afz8/IIXLVqUqRri4uLInz9/5hqQizliux2xzWC7dnvv3UuSry+3nnySw4ej+eakE9dc8tKpohsNn3TBSRk3TIZsa+s0aNAgXGtdM90ZtdbZ+gB8gc1AYcAVWAm8Bhy/Z54SwP701hUcHKwzKzQ0NNPL5maO2G5HbLPWNmr3unVau7vrpIaN9fAle3TJEWt1y0//0MeiY7P+vTJBtrV1gDBtxee1EbuYGgOntNaXAZRSK4BnAB+llIvWOgXwB84bUJsQ4mHWrUO3acOVkuVoV6M3UX9HMbBBOYY0Li93d7NTRgTEGaCuUiovpl1MjYAwIBRoj+lMpu7AKgNqE0JYoA8eJKVDR04UKknHF8dQO6AsXzatQKWicl2DPcv2gNBa71JKLcN0KmsK8A/wFfAzsEgp9aF52tzsrk0I8aCoawlE9RhCWVyYNmgy87vWp/qTvkaXJbKBIWcxaa3HAGPum3wSqG1AOUIIC1JS0/j+z9NM/vUIbo0H8f4beZndubHcq8GByI5DIcQDwk9fo/0nm4kfPoJn/NxY/VYTWndpIuHgYGSoDSHEXRdvJDJlwxG2he7h21UfUen8UfSwDqgCeY0uTRhAAkIIQUJSCl/+fpKvtp6keuQ+fls7iXypSbBqFap1a6PLEwaRgBDCwW08FM2YVQc4fyORD2N20/WHD1Bly5ruHV25stHlCQNJQAjhoC7cuMXY1QfZsucMtfKlMuP1RtRyDoR8N+GDD8ABr0gW/yYBIYSDSUvT/LjrNHunfEHXvb/x+dkDONV7HqdxHYACMG2a0SWKHEICQggHcvJyHFO+3kCnOeN59dTfJJcug8vr/aBlS6NLEzmQBIQQDiD+dgozQ48z949TdN+zmWcuHkZ/+imuAwaAk5ztLiyTgBDCzq3Ze54Jq/eTP/I4LV94lj5vT8I1diT4+xtdmsjhJCCEsFO3UzVvL9vLuu1HmLd+KtXOReD80THwzmt6CJEOCQgh7NDxS3F8sPMWN68eZtPqsRQ+cxw1ezb4+RldmshFJCCEsCNpaZrvd51mwi+H8UhJYuuOz/A5ewLWroUXXjC6PJHLSEAIYSfOxiTw9rJ97Dx5lfoVCvPmrnn47N4O338v4SAyxaqAUEpVwHR/aD+tdYBSKhBorbX+0KbVCSHSlZicytxtp/g89DhOSvFxu6p0rFmCrSVj4eXW0KyZ0SWKXMra89u+BkYByQBa631AZ1sVJYRIn9aa9Qcu0mTa70z+9QjPlSvE5ufc6fR2d9TFi2hnZwkH8Vis3cWUV2u9W/37ZuQpNqhHCGGFA+du8MHaQ+w6FUMFv/z82LkKz8yeAF99BcWKwalTRpco7IC1AXFFKVUW0ABKqfbABZtVJYSwKCY+ifE/R7Dinyh887rxQZsqvJIvFpdOzeHoUXjzTRgzBjw9YcsWo8sVuZy1ATEA021Bn1JKnQNOAa/arCohxANCj1zi7WX7uJ6QRN96ZRjQoBxe7q7Qowdcuwa//QYNGhhdprAjVgWE1vok0FgplQ9w0lrH2rYsIcQdCUkpTPjlMN/9eZqKfp4s+L/aVPLLbwoF94Lw2WcQFwdFixpdqrAzVh2kVkp9pJTy0VrHa61jlVK+Sik5g0kIG9tx/AovTN/Kd3+epvdzpVk18Fkq+bpBly4QEgIJCabdSRIOwgasPYupudb6+p0nWutrQAvblCSEiE1MZtSK/XSZswtnpVjS72n+27Iy7uejoGlTWLwYuncHDw+jSxV2zNpjEM5KqTxa69sASikPII/tyhLCcYVFxvDGkj2cu3aLvvXK8EbjCnioNPjwQ/joI9NMCxdCZznTXNiWtQHxPbBJKfUtpjOZ/g+Yb7OqhHBAyalpfLbpGDNDj1Pc14Ml/Z6mZqkCphdTMd0C9MUXYepUePJJQ2sVjsHag9STlFL7gUaAAj7QWv9q08qEcCCRV+IZsngPe89ep22N4rzfugqerk4wZQr07AkFC0JoqNwGVGQrq8di0lqvA9bZsBYhHI7WmqVhUYxdcxAXJ8XMLtVpGVgMLl407UL6/XfIkwcGDZJwENnO2rGY2gIfA09g6kEoQGutvWxYmxB27dLNREavPMDGQ9HULVOATzpWo5iPBxw5Yhoi49IlmD8funUzulThoKztQUwCWmmtI2xZjBCOQGvNyj3nGLv6EInJqbzT4il6PVcGZycFYWGmkVddXEy9h5o1jS5XODBrAyJawkGIxxd9M5HRP+3nt4hLBJf0ZVL7QMoWzg9paYAy3Qa0Zk2YPRvKlDG6XOHgrA2IMKXUYmAlcPvORK31CptUJYSd0Vqz/O9zjFtzEJ2YyLTyTrQulYjzto2wdavpsXs3FCkCv8r5HyJnsDYgvIAEoOk90zQgASHEw8THQ0wMJ1NcGLPpNH8cv0qtUr58ffg3fHpP+d98Li6mA9KxseDjY1y9QtzH2tNce2blmyqlfIA5QAD/u67iCLAYKAVEAh3NV2wLkbtcvAgjR6KXLkUlJFAGmI9i0yfzaNS3BU5ni8MztaBwYfD1Ne1WKlTI6KqFeIC1ZzG5A72AKoD7nela6//L5PvOANZrrdsrpdyAvMA7wCat9USl1EhgJDAik+sXIvvdvm06JdXTk9vrN7C+SgN2FijNs4VdaVjEjSaVngAnBSVLmh5C5HDW7mL6DjgMvACMA7oCmTporZTyAuoBPQC01klAklKqDRBinm0+sAUJCJEbxMfDe+/Bhg3c3LaTCZtOsrjbF5T28+LjdoH/uxpaiFzG2oAop7XuoJRqo7Wer5T6EcjskbQywGXgW6VUEBAODMF0v+sLAFrrC0qpJzK5fiGyz7Fj0LIlHD3K8bav8n9TQ4lKdqJvgwoMbVwed1dnoysUItOU1jr9mZTarbWurZTaCvQHLgK7tdYZPg9PKVUT+BN4Vmu9Syk1A7gJDNJa+9wz3zWtta+F5fsCfQH8/PyCFy1alNESAIiLiyO/A16Z6ojttlWbvffsIWDMGFJQjGg/ktWFq1LOx4kuldwo4218MMi2dhwZbXeDBg3CtdbpX2SjtU73AfQGfDHtGjoJXAL6WbOshXUVASLvef488DOmg9RFzdOKAkfSW1dwcLDOrNDQ0Ewvm5s5Yrtt0ua0NH27dh19vlhp/Vy/Ofq5jzfpNXvP6bS0tKx/r0ySbe04MtpuIExb8Xlt7S6mTdp0RtFWTLuIUEqVtnLZ+wPpolLqrFKqotb6CKYBAA+ZH92BieafqzKzfiFsKimJ1FuJzNlzmR+eGUyCsyvdm1ejT70ysjtJ2B1rA2I5UOO+acuA4Ey+7yDgB/MZTCeBnphuXrREKdULOAN0yOS6hbCNCxe4/XI79iU4MaH5SBpXq8CYVpUpUSCv0ZUJYROPDAil1FOYTm31Ng/Yd4cX95zumlFa6z2Apf1fjTK7TiFsJi0N/e233H57JGmxsSxuPYxpnavxUrXiKKWMrk4Im0mvB1ERaAn4AK3umR4L9LFVUULkGEePktihE+779nCoWEV+HDCNoW+0xd9Xeg3C/j0yILTWq4BVSqmntdY7s6kmIXKExORU5uy/TsPomyxoO4Kqb73OpNolcXKSXoNwDNYeg3hZKXUQuAWsB4KAoVrr721WmRBG0ZpDM+YyMLEkJ68ncXzaUka3rEJhT7kNu3AsTlbO11RrfRPT7qYooALwls2qEsIgp6JvsrFNTyq/0YdWu3/hh951mP5KDQkH4ZCs7UG4mn+2ABZqrWPk4JywJwfO3eC3WYtpOmciTS6dYm+rV/jPwkm453FNf2Eh7JS1AbFGKXUY0y6m/kqpwkCi7coSInvE3U5hwi8ReE+fwttbF3DjiaLc/PY7grp3BfkSJByctcN9j1RKfQzc1FqnKqXigTa2LU0I29q3ajMfbj/PX86+jOzcjsSQ0niPHAEeHkaXJkSOkN51EA211pvvvQbivl1LcsMgketcik0kvO9wmi/6nG61W/D24u/MI662Nro0IXKU9HoQ9YDNmK6B0IC676cEhMg1bqekMvePkzi/9x79ti3iQEhLGi6ZR97CMhy3EJakFxCxSqlhwAH+FwyYfxci1zh4/gZvLt5Du4XT6fPXSmJf7UHA/LngZO2JfEI4nvQC4s74sRWBWpgG0FOYehRbbViXEFkiJTWN1cdvs3rDdgp4uPBKzCEYNAjP6dMlHIRIR3pXUr8PoJTaANTQWsean48Fltq8OiEyKjUVpkyBpUtJTNVE3khiwI3rqCmL+G/HWuR/5WcoW1bOUBLCCtZ+hXoSSLrneRJQKsurESIjUlIgNvbf01avhpEjuZiYxu6biuvKleTyZZja0B/ffG5QrpyEgxBWysg9qXcrpX7CdPzhZUz3jRbCGBs2QL9+EBkJZcpAjRqwaBERdRqyeMAnzMtXniZVivDRy1U5GL6Tp0qWNLpiIXIda6+DGK+UWofp7m8APbXW/9iuLCEeYflyaN8eKlaEsWPhwAFSr8YwecMx5vxxEi+/AKa3rEybasVkOG4hHoO1PQi01n8Df9uwFiGs8+KLMHkyDBwI7u5sP36FUSv2c+b3E3Ss6c+o5pVMu5OEEI9FTuMQuUNKCgwbBlFR4O4Ow4dzPc2Jt5bupeucXTg7KRb2qcuk9kESDkJkEat7EEIYavhwmDEDgoKge3fCT8cw8Md/uBR7m/4hZRncqLzcE1qILCYBIXK+uXNN4TBkCLpbN77eeoJJ649QzMeDlf2fpaq/t9EVCmGXJCBEzrZ+Pbz+OjRpwqUxHzJqfhibDl+ieUARPm4fiJe7DMcthK1IQIicS2t4/30ICGDj+zN5+7MdJCSlMrZVZbo/U0rOUBLCxiQgRM6lFDeWrWTKmn18t+oYgf7efNIxiHJPeBpdmRAOQQJC5DznzsGECWzo+SbvrDvB9YRkhjQqz8CG5XB1lhPvhMguEhAiZ4mJIaVJU1JORTLlViX8ggJZ8H+1qVzMy+jKhHA4EhAix0iNi+dq/Sb4HD1G387v07prU/rVLyu9BiEMIgEhjLNzJ7i5QXAwpw9HQsOGlLgQyWf/Gc+4cYMpVSif0RUK4dDkq5kwxpo10KgResgQlv51hn4zNxOPM3vGTWPw5yMlHITIAaQHIbLfnDnQrx8p1Wvw7qvvs3D5fuoGVsJ3/D6KensYXZ0QwkwCQmSf+HjTKKzr1xP9TAjtGrzBxTO3Gd60Av8JKYezk1zXIEROIgEhsk/evCR45GN1x8H898mGVClakDntqvJUETlDSYicyLCAUEo5A2HAOa11S6VUaWARUADTsOKvaa2THrUOkUuMGUNq23Z8F+fJ5IDeaGBU04r0eKaU9BqEyMGM7EEMASKAO18fPwamaa0XKaW+AHoBs40qTmSRuXNh3DiWhEUxtmpb6lcozIcvBVCiQF6jKxNCpMOQs5iUUv7Ai8Ac83MFNASWmWeZD7xkRG0i6+h//iG1f3+2l67G5Frtmd6pGvN61pJwECKXMKoHMR14G7gzqE5B4LrWOsX8PAoobkRhImtcPH0e1aw1aW6eLBwykfU96/OEl7vRZQkhMkBprbP3DZVqCbTQWvdXSoUAw4GewE6tdTnzPCWAX7TWVS0s3xfoC+Dn5xe8aNGiTNURFxdH/vz5M9eIXMzW7U7TmtCzKTw5+0u6717F/FGTKdcoyNCRV2VbOw5HbDNkvN0NGjQI11rXTHdGrXW2PoAJmHoIkcBFIAH4AbgCuJjneRr4Nb11BQcH68wKDQ3N9LK5mS3bfeTiTd121nZdcsRa3WP2Vn3hl002e6+MkG3tOByxzVpnvN1AmLbi8zrbdzFprUcBowDu9CC01l2VUkuB9pjOZOoOrMru2kTmxN9O4dNNx/hxcwRv7fiBHh+MpWXDALlfgxC5XE66DmIEsEgp9SHwDzDX4HpEOtLSNGv2nee7Hzbzwu8r2HVoE3njYyGmM6gH9g4KIXIZQwNCa70F2GL+/SRQ28h6hHW01mw5cpnJ6yLosuBjluxdD87OOLVrB0OHQt26RpcohMgCOakHIXKB2MRk3lq6j/UHL/Jkgbw8X8obVX8w6q23oLiceCaEPZGAEFY7Gh3L69+F02r1HBp268RLPerj5hwCcqxBCLskASHSpbVm+d/neHflATod2swbf/wAIWXApZXRpQkhbEgCQjzSlbjbvLNiPxsORdM95TRj1kyHhg3h3XeNLk0IYWMSEMIirTW/7L/Ie6sOEJuYwsRAdzoNehdVpgwsWwaurkaXKISwMQkI8YDTV+N5b9VBfj96mYDiXizsWI0K7w03vbh2Lfj6GlugECJbSECIu5JS0vhq6wk+23wcV2cn3mtZmW7BRXHxcIdp02DYMChXzugyhRDZRAJCALD37HV+mDiPOlvXsDotllJ+XuTZoeDyZdi6Fby8oGJFo8sUQmQjCQhHdv48yYsW82nxuny+J4aBp4/RKmoPbpWfglvxkJpquuhNjjcI4ZAkIBzR7t3w3nvoDRtw1ZrjbUbSuXsXeo+YiJvXZ+BkyG1ChBA5jASEI7l1iyrvvQd//EG8ly9fP9OZXbWbMrB/S54tV8jo6oQQOYwEhCPx8CAhWTP/hZ5MqtScNvWe4usWlcifR/4ZCCEeJJ8M9u72bXjnHVL7D+CLM2l8Un8EBT3d+bx9ICEVnzC6OiFEDiYBYc9iY6FZM9ixg2/PKyaXbEDtIi582acevvncjK5OCJHDSUDYq+RkaN+etF27eLv9O6wr/zygKlbWAAASQElEQVRT2wRQ4OYxCQchhFXkdBV7pDVxPXrBhg2MbDKAYyEt+GXI87QL9pe7vAkhrCY9CDuTlJLGN+v3Ufv3MHY9/wrlRw1h/LOlcHWW7wJCiIyRgLAX8fFcmPkV4694s9bZjzYffsPIlwIp6pPX6MqEELmUBERut28fqV99TfL8BRSNu0mNeq/QetYnNK1SxOjKhBC5nAREbnXrFjRpAtu3k+rsyoYKT3O8Q3d6vNWVAvnzGF2dEMIOSEDkUifjUvmtygtccK/M38+/yJBOdRn2lJ/RZQkh7IgERC4Tu3Eza7dG8G5qadyL1GZQly4sfbY0bi5yEFoIkbUkIHIJrTU7P/6S6u8OJahAMTp9voKhL1SisKfsThJC2IYERC5w/kos/7zWnxfXf8/R0lVwXb2a8QFljC5LCGHnJCBysLQ0zY+7TlO0Z1dePLKDI21fpfwPc3Byl16DyN2Sk5OJiooiMTExS9fr7e1NRERElq4zN3hYu93d3fH398c1k/d0kYDISWJjYd8+2LOH63sP0qfGq/wVeY1vCxXg6qAvqTigr9EVCpEloqKi8PT0pFSpUll6dX9sbCyenp5Ztr7cwlK7tdZcvXqVqKgoSpcunan1SkDkBFrDm2+a7vtslpzPl0s+DZjcsS4h41einJ0NLFCIrJWYmJjl4SD+TSlFwYIFuXz5cqbXIQGRE2zdCtOmcaNjF2b4BPKLcxGCnw1gWesAOQgt7JaEg+097t9Yzo3MAVKfr8eaad9Rq2xXVvrX4N3Xm/J512AJByFyiHnz5jFw4MBMLdujRw+WLVuW7vrPnz9/93nv3r05dOhQpt4vK0kPwihXrsDMmVys8zyDTnvw10VfXqjyBONfrkohuRJaCIcyb948AgICKFasGABz5swxuCKTbO9BKKVKKKVClVIRSqmDSqkh5ukFlFIblVLHzD99s7u2bHHyJPTrhy5RAt5/n23vTObwxVimdQrii1eDJRyEyEYLFiwgMDCQoKAgXnvtNdasWUOdOnWoXr06jRs3Jjo6+oFloqOjefnllwkKCiIoKIgdO3YQGRlJQEDA3XmmTJnC2LFjH1h23Lhx1KpVi4CAAPr27YvWmmXLlhEWFkbXrl2pVq0at27dIiQkhLCwMAAWLlxI1apVCQgIYMSIEXfXlT9/fkaPHk1QUBANGza0WOvjMqIHkQK8qbX+WynlCYQrpTYCPYBNWuuJSqmRwEhgxCPWk/ts2EBax46k3kpkVUBDvqjeCv9ng/n15aoU8/EwujohDPH+moMcOn8zS9aVmpqKs7MzlYt5MaZVlUfOe/DgQcaPH8/27dspVKgQMTExKKX4888/UUoxZ84cJk2axNSpU/+13ODBg6lfvz4//fQTqampxMXFce3aNavqGzhwIO+99x4Ar732GmvXrqV9+/bMnDmTKVOmULNmzX/Nf/78eUaMGEF4eDi+vr40bdqUlStX8tJLLxEfH0/dunUZP348Q4cO5euvv+a///1vBv5a6cv2gNBaXwAumH+PVUpFAMWBNkCIebb5wBbsKCDS0jQ7dh+hkHtBenf5L0HPBvFJ/TIE+vsYXZoQDmnz5s20b9+eQoUKAVCgQAH2799Pp06duHDhAklJSRZPD928eTMLFiwAwNnZGW9vb6sDIjQ0lEmTJpGQkEBMTAxVqlShVatWD53/r7/+IiQkhMKFCwPQtWtXtm7dyksvvYSbmxstW7YEoFq1amzbti1D7beGoccglFKlgOrALsDPHB5orS8opZ4wsLSsEx7OxW1/MdAtkLC4Mjw3Zj7z2lWj3BP5ja5MiBwhvW/6GZGR6yC01g+c5TNo0CCGDRtG69at2bJli8XdRJa4uLiQlpZ297mlCwATExPp378/YWFhlChRgrFjx6Z7oaDW+qGvubq63q3f2dmZlJQUq2rNCMMCQimVH1gODNVa37T2dCylVF+gL4Cfnx9btmzJ1PvHxcVlellrFNq2jSe//x6vI0dw8/Di+JD59ArIx3PFk4k6FEaUQSco2LrdOZEjthlydru9vb2JjY3N8vWmpqZavd66devSpUsXevfuTcGCBYmJieHatWv4+PgQGxvLnDlz7q4vMTGRpKQkYmNjqVevHtOmTWPAgAGkpqYSHx9P3rx5iY6OJjIykvz587Nq1SoaN25MbGwsycnJ3Lp1i8uXL6O1Jk+ePFy4cIElS5bQpk0bYmNj8fDwIDo6+m7td9ZbpUoVBg8eTGRkJD4+Pnz//ff069fv7nx3fqalpZGcnGyx7YmJiZn+d2BIQCilXDGFww9a6xXmydFKqaLm3kNR4JKlZbXWXwFfAdSsWVOHhIRkqoYtW7aQ2WXTtXUreswYogqXYErjfsS178xvXZ/OEQegbdruHMoR2ww5u90RERE2ueI5Iz2I2rVr8+6779KyZUucnZ2pXr0648aNo0ePHhQvXpy6deveveLb3d0dNzc3PD09mTVrFn379uWHH37A2dmZ2bNn8/TTTzNmzBgaN25M6dKlqVKlCnny5MHT0xNXV1c8PDwoUaIEffv25ZlnnqFUqVLUqVPn7jy9e/dm2LBheHh4sHPnTpydncmXLx/ly5dn4sSJtGrVCq01LVq0oHPnznfbcKetTk5OuLq6Wmy7u7s71atXz9wfVGudrQ9AAQuA6fdNnwyMNP8+EpiU3rqCg4N1ZoWGhmZ62Ue5cCVWX/Hz1yd9i+lnR6/Ua/aes8n7ZJat2p2TOWKbtc7Z7T506JBN1nvz5k2brDene1S7Lf2tgTBtxee1ET2IZ4HXgP1KqT3mae8AE4ElSqlewBmggwG1ZVr87RQ+Dz3ON9tP8VSL4bwQ5M+qXs0omAN6DUIIkRlGnMW0DVMvwpJG2VlLVtBa8/P+C3y4NgKPyBM0a1iLN9+oT4kCeY0uTQghHotcSf0YzsYkMGrFfrYdv0Kb5PNM/3YAqu6XUCCT+/uEECIHkYDIhLQ0zfe7TjNx3WGclGJ8s7J0+c9wlJ8fvPyy0eUJIUSWkIDIoMgr8YxcsY8/T8ZQr0JhJrStSvExI+HwYdi4EXztc4QQIYTjkYCwUmqa5pttp5i68QiuTk583K4qHWuWQG3aBJ9+CoMGQePGRpcphBBZRob7tsLhizdpO3sH43+J4Llyhdg4rD6daj1puorx+nWoXRsmTjS6TCFEBn366adUqlQJX19fJpr/D48dO5YpU6YADw7D7WikB/EIicmpfLb5GF/+fhIvD1dmdK5G66Bi/748v317aNcO5OYnQuQ6s2bNYt26dQ+9Jef9w3BbIyUlBRcX+/holR7EQ+w9e50WM/7g89ATtKlWnE3D6tOmWvH/hUNoKHzxBaSlSTgIkQu9/vrrnDx5ktatWzNt2rQHbghkaRju8PBw6tevT3BwMC+88AIXLlwAICQkhHfeeYf69eszY8YMI5pjE/YRc1koLU0zZ9tJJq0/gp+XO9/1qs3zT3pBHjfTDImJkJICvXqBiwt07w4eMlS3EI/F0pAgHTtC//6QkAAtWjz4eo8epseVK6aePOCRmgrOzmDF2ENffPEF69evJzQ0lLVr1z7w+v3DcCcnJzNo0CBWrVpF4cKFWbx4MaNHj+abb74B4Pr16/z+++/WtzkXkIC4x+XY2wxfuhfvlcto0+wF3utaF++P3ofNm2HbNjh4EJo1g4AAiIw03UtawkEIh3DkyBEOHDhAkyZNANOAekWLFr37eqdOnYwqzWYkIK5dg/HjiTpxjtY1/g+f6Cjmr5mMXj8d9cmTcOIE9OsHWoO7OwQGwqZNprOWnnvO6OqFsA+P+safN++jXy9U6O7rtzIwWF9Gaa2pUqUKO3futPh6vnz5bPK+RnLsYxBHj5Jaqzap06az+1g0xb3c+PLd9rBvH2rgQChQABYuNB1rcHeHypVN1zocOACffGJ09UIIG/P09Lw7hHbFihW5fPny3YBITk7m4MGDRpZncw7bg9A7wkmY0IbENOjXZQJ1Xm3N8kblcXNxAr+qDw8ApaBK1t3gRAiRc/Xo0YPXX3/97jDcy5YtY/Dgwdy4cYOUlBSGDh1KFTv+PFD6EXcsyulq1qyp79zYOyO+3XiQF19+nuseXqwYO4v2HZ6n3BO26ZbmNDn5HgG24ohthpzd7oiICCpVqpTl683I/SDsyaPabelvrZQK11rXtLjAPRyyBxFcqTjfDP2QXv9pz8ji9nFnUyGEyGoOGRCB/j7ENK5MYQkHIYR4KMc+SC2EEOKhJCCEEIbIzcc/c4vH/RtLQAghsp27uztXr16VkLAhrTVXr17F3d090+twyGMQQghj+fv7ExUVxeXLl7N0vYmJiY/1gZhbPazd7u7u+Pv7Z3q9EhBCiGzn6ur60BFUH8eWLVuoXt3xbvlrq3bLLiYhhBAWSUAIIYSwSAJCCCGERbl6qA2l1GXgdCYXLwRcycJycgtHbLcjthkcs92O2GbIeLtLaq0LpzdTrg6Ix6GUCrNmLBJ744jtdsQ2g2O22xHbDLZrt+xiEkIIYZEEhBBCCIscOSC+MroAgzhiux2xzeCY7XbENoON2u2wxyCEEEI8miP3IIQQQjyCQwaEUqqZUuqIUuq4Umqk0fXYglKqhFIqVCkVoZQ6qJQaYp5eQCm1USl1zPzT1+has5pSylkp9Y9Saq35eWml1C5zmxcrpdyMrjGrKaV8lFLLlFKHzdv8aQfZ1m+Y/30fUEotVEq529v2Vkp9o5S6pJQ6cM80i9tWmXxq/mzbp5Sq8Tjv7XABoZRyBj4HmgOVgVeUUpWNrcomUoA3tdaVgLrAAHM7RwKbtNblgU3m5/ZmCBBxz/OPgWnmNl8DehlSlW3NANZrrZ8CgjC13663tVKqODAYqKm1DgCcgc7Y3/aeBzS7b9rDtm1zoLz50ReY/Thv7HABAdQGjmutT2qtk4BFQBuDa8pyWusLWuu/zb/HYvrAKI6prfPNs80HXjKmQttQSvkDLwJzzM8V0BBYZp7FHtvsBdQD5gJorZO01tex821t5gJ4KKVcgLzABexse2uttwIx901+2LZtAyzQJn8CPkqpopl9b0cMiOLA2XueR5mn2S2lVCmgOrAL8NNaXwBTiAD2dt/V6cDbQJr5eUHgutY6xfzcHrd3GeAy8K1519ocpVQ+7Hxba63PAVOAM5iC4QYQjv1vb3j4ts3SzzdHDAhlYZrdnsqllMoPLAeGaq1vGl2PLSmlWgKXtNbh9062MKu9bW8XoAYwW2tdHYjHznYnWWLe794GKA0UA/Jh2sVyP3vb3o+Spf/eHTEgooAS9zz3B84bVItNKaVcMYXDD1rrFebJ0Xe6nOafl4yqzwaeBVorpSIx7TpsiKlH4WPeBQH2ub2jgCit9S7z82WYAsOetzVAY+CU1vqy1joZWAE8g/1vb3j4ts3SzzdHDIi/gPLmMx3cMB3UWm1wTVnOvO99LhChtf7knpdWA93Nv3cHVmV3bbaitR6ltfbXWpfCtF03a627AqFAe/NsdtVmAK31ReCsUqqieVIj4BB2vK3NzgB1lVJ5zf/e77Tbrre32cO27Wqgm/lsprrAjTu7ojLDIS+UU0q1wPTN0hn4Rms93uCSspxS6jngD2A//9sf/w6m4xBLgCcx/QfroLW+/wBYrqeUCgGGa61bKqXKYOpRFAD+AV7VWt82sr6sppSqhunAvBtwEuiJ6QugXW9rpdT7QCdMZ+39A/TGtM/dbra3UmohEIJpxNZoYAywEgvb1hyUMzGd9ZQA9NRah2X6vR0xIIQQQqTPEXcxCSGEsIIEhBBCCIskIIQQQlgkASGEEMIiCQghhBAWSUAIYSXziKn9zb8XU0otS28ZIXIzOc1VCCuZx7Raax45VAi755L+LEIIs4lAWaXUHuAYUElrHaCU6oFpNE1nIACYiumCtdeA20AL80VMZTENNV8Y00VMfbTWh7O/GUJYR3YxCWG9kcAJrXU14K37XgsAumAaTn48kGAeOG8n0M08z1fAIK11MDAcmJUtVQuRSdKDECJrhJrvuxGrlLoBrDFP3w8EmkfVfQZYahoNAYA82V+mENaTgBAia9w71k/aPc/TMP0/c8J0n4Jq2V2YEJklu5iEsF4s4JmZBc334jillOoAd+8dHJSVxQmR1SQghLCS1voqsN188/jJmVhFV6CXUmovcBA7vNWtsC9ymqsQQgiLpAchhBDCIgkIIYQQFklACCGEsEgCQgghhEUSEEIIISySgBBCCGGRBIQQQgiLJCCEEEJY9P8dy38cYVfzAAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1087c4588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGcxJREFUeJzt3Xt0XeV55/Hvo6styZZsS75LlnGwwAR8QSYGQm0gwy0X2g5pSzKQZGB5WGVamGFm0mamyepK+0dWUpJ2UsJyIXXaMqRTcANlcmnGmFDG4Fo2xjfZxuC7DZZvsi3fdHnmj7NtC/kcnSNpH2/tfX6ftbSsc857zn422/zO63e/e7/m7oiISLIURV2AiIiET+EuIpJACncRkQRSuIuIJJDCXUQkgRTuIiIJpHAXEUkghbuISAIp3EVEEqgkqg3X1tZ6Y2NjVJsXEYmlNWvWHHL3umztIgv3xsZGWlpaotq8iEgsmdmuXNppWEZEJIEU7iIiCaRwFxFJIIW7iEgCKdxFRBJI4S4ikkAKdxGRBIpduG/94ATf+cVWDp88G3UpIiLDVuzC/f22k3x/xXYOnlC4i4hkErtwryhPXVR76lxXxJWIiAxfsQv3qvJiAE6e7Y64EhGR4St24V5RFvTcz6rnLiKSSezCvTII945z6rmLiGQSu3CvCIZlNOYuIpJZ7ML9Qs9dY+4iIhnFLtxHlBZRZNChMXcRkYxiF+5mRmVZCR0alhERySh24Q6pcfdTGpYREckoa7ibWb2ZrTCzVjPbZGaPZWi3yMzWBW1+FX6pF6nnLiLSv1zWUO0CnnD3tWY2ClhjZr90983nG5hZDfAUcJe77zaz8XmqFwh67poKKSKSUdaeu7sfcPe1we8ngFZgSp9mXwCWufvuoN3BsAvtraKsRCdURUT6MaAxdzNrBOYCq/q8NBMYY2avmdkaM3sww/sXm1mLmbW0tbUNpl4Aqso1LCMi0p+cw93MqoAXgcfd/Xifl0uA64FPA3cCf2RmM/t+hrsvcfdmd2+uq6sbdNEVZTqhKiLSn1zG3DGzUlLB/py7L0vTZC9wyN07gA4zex2YDWwLrdJedEJVRKR/ucyWMeBZoNXdn8zQ7CXgFjMrMbMK4BOkxubzQlMhRUT6l0vP/WbgAWCDma0Lnvsa0ADg7k+7e6uZ/RxYD/QAz7j7xnwUDBd77u5O6rtHRER6yxru7v4GkDVB3f3bwLfDKCqbivJiehzOdvUworT4cmxSRCRWYnmFalWwGtNJTYcUEUkrluF+ccEOjbuLiKQTy3CvLEsNxWjGjIhIerEMdy2SLSLSv1iG+4Weu4ZlRETSime4l59fjUk9dxGRdOIZ7lokW0SkX7EMdy2SLSLSv1iGuxbJFhHpXyzD/fwi2eq5i4ikF8twv7BItnruIiJpxTLcITXurtkyIiLpxTbcdU93EZHMYhvuWiRbRCSz+Ia7FskWEckotuFeWaaeu4hIJvEN93L13EVEMolvuOuEqohIRrkskF1vZivMrNXMNpnZY/20nW9m3WZ2X7hlXkqLZIuIZJbLAtldwBPuvtbMRgFrzOyX7r65dyMzKwa+BfwiD3VeQotki4hklrXn7u4H3H1t8PsJoBWYkqbp7wEvAgdDrTCD3otki4jIRw1ozN3MGoG5wKo+z08BfgN4OqzCsqnSPd1FRDLKOdzNrIpUz/xxdz/e5+XvAV91934Hwc1ssZm1mFlLW1vbwKvtpUJ3hhQRySiXMXfMrJRUsD/n7svSNGkGfhyMfdcC95hZl7v/pHcjd18CLAFobm72oRSuRbJFRDLLGu6WSuxngVZ3fzJdG3ef3qv9UuCVvsEeNi2SLSKSWS4995uBB4ANZrYueO5rQAOAu1+2cfbetEi2iEhmWcPd3d8Acp5r6O5fHkpBuTo/5q6eu4jIpWJ7herF2TLquYuI9BXbcD+/SLZOqIqIXCq24a5FskVEMottuI8oLcK0SLaISFqxDXctki0ikllswx2gsrxYPXcRkTTiHe5lJZzUvWVERC4R63DXItkiIunFO9y1SLaISFqxDnctki0ikl6sw72iXOuoioikE+twryor0TqqIiJpxDrcK8qLNeYuIpJGrMO99yLZIiJyUazDXYtki4ikF+twv3jzMA3NiIj0Fu9wv7DUnk6qioj0Fu9wD5ba0y0IREQ+Ktbhfn6RbIW7iMhHZQ13M6s3sxVm1mpmm8zssTRtvmhm64OflWY2Oz/lflTjuAoAtn144nJsTkQkNnLpuXcBT7j71cAC4FEzm9WnzQ5gobtfB3wTWBJumek1jK1gbGUZ63YfuxybExGJjZJsDdz9AHAg+P2EmbUCU4DNvdqs7PWWt4CpIdeZlpkxe2o16/Yo3EVEehvQmLuZNQJzgVX9NHsI+NngSxqYOfVj2N52khNnOi/XJkVEhr2cw93MqoAXgcfd/XiGNreSCvevZnh9sZm1mFlLW1vbYOq9xJyGGtxh/d72UD5PRCQJcgp3MyslFezPufuyDG2uA54B7nX3w+nauPsSd2929+a6urrB1vwRs6dWA2hoRkSkl1xmyxjwLNDq7k9maNMALAMecPdt4ZbYv5qKMqbXVircRUR6yXpCFbgZeADYYGbrgue+BjQAuPvTwNeBccBTqe8Cuty9Ofxy05tTX8Mb2w/h7gTbFxEpaLnMlnkD6Dcx3f1h4OGwihqoOfU1/OPb+9jffoYpNSOjKkNEZNiI9RWq582prwHQfHcRkUAiwv3qSaMpKyli3Z6jUZciIjIsJCLcy0qKuGbyaJ1UFREJJCLcITU0s2FfO53dWrhDRCRR4X6ms4etH+gmYiIiiQn3G6aPxQyWtx6MuhQRkcglJtwnVY/kphnjeGHtHnp6tGC2iBS2xIQ7wOevr2fPkdOs2nEk6lJERCKVqHC/85qJjCov4R/W7Im6FBGRSCUq3EeWFfOZ2ZP42YYPtPSeiBS0RIU7wH3X13O6s5ufrj8QdSkiIpFJXLjPa6jhirpKDc2ISEFLXLibGfddP5XVO4+y41BH1OWIiEQiceEO8JtzU0u4/p/1+yOuREQkGokM94nVI/j4lNG8tjWcpfxEROImkeEOcGvTeNbuPkr7KS2cLSKFJ7Hhvqipjh6Hf9mu3ruIFJ7Ehvuc+jFUjyxlxRaFu4gUnsSGe3GRccuVtfxqW5vuNSMiBSdruJtZvZmtMLNWM9tkZo+laWNm9hdmtt3M1pvZvPyUOzC3No3n0MmzbD5wPOpSREQuq1x67l3AE+5+NbAAeNTMZvVpczdwZfCzGPhBqFUO0q/NrAPgta26DbCIFJas4e7uB9x9bfD7CaAVmNKn2b3A33jKW0CNmU0KvdoBqhtVzrVTqlmhKZEiUmAGNOZuZo3AXGBVn5emAL2v99/LpV8AmNliM2sxs5a2tssTuLc21fH27qMcO3XusmxPRGQ4yDnczawKeBF43N37DmJbmrdcchbT3Ze4e7O7N9fV1Q2s0kFa2DQ+NSXy3UOXZXsiIsNBTuFuZqWkgv05d1+WpsleoL7X46nAsLj2f059DaPKS3jz/cNRlyIictnkMlvGgGeBVnd/MkOzl4EHg1kzC4B2dx8W99wtLjLmTRvDaq3OJCIFJJee+83AA8BtZrYu+LnHzB4xs0eCNj8F3ge2A38F/G5+yh2cG6aP5d2DJznaoXF3ESkMJdkauPsbpB9T793GgUfDKipszdPGANCy6yj/ZtaEiKsREcm/xF6h2tvs+hrKioto2amhGREpDAUR7iNKi7l2ajWrFe4iUiAKItwB5jeOZcO+ds50dkddiohI3hVQuI+hs9tZt+dY1KWIiORdwYR787SxAJoSKSIFoWDCvbqilKYJo1i962jUpYiI5F3BhDvA/OljWLvrKN26v7uIJFxhhXvjWE6e7aJV93cXkYQruHAHNCVSRBKvoMJ9cs1IJlePYI3G3UUk4Qoq3AHmTUuNu4uIJFnhhXvDGPa3n+FA++moSxERyZuCC/frg5uIrd2li5lEJLkKLtxnTR7NiNIijbuLSKIVXLiXFhdx3ZQa1u5WuItIchVcuEPqpOqm/bqJmIgkV2GGe0MNnd3Oxn3tUZciIpIXhRnuwUlVjbuLSFLlskD2D83soJltzPB6tZn9k5m9Y2abzOwr4ZcZrtqqchrHVWjcXUQSK5ee+1Lgrn5efxTY7O6zgUXAn5lZ2dBLy695DWNYs+sYqeVfRUSSJWu4u/vrQH83Y3FglJkZUBW07QqnvPyZN20Mh06eZc8RXcwkIskTxpj794Grgf3ABuAxd+8J4XPzal5DcDGThmZEJIHCCPc7gXXAZGAO8H0zG52uoZktNrMWM2tpa2sLYdOD1zRxFFXlJbTs0h0iRSR5wgj3rwDLPGU7sAO4Kl1Dd1/i7s3u3lxXVxfCpgevuMiY21BDy0713EUkecII993A7QBmNgFoAt4P4XPzrnnaWLZ+eIL2051RlyIiEqpcpkI+D7wJNJnZXjN7yMweMbNHgibfBG4ysw3AcuCr7n4ofyWHZ37jGNw17i4iyVOSrYG735/l9f3AHaFVdBnNaaihpMho2XmEW5vGR12OiEhoCvIK1fMqykq4Zko1qzXuLiIJU9DhDjB/2hje2XOMs126iZiIJEfBh3tz41jOdvWwcd/xqEsREQmNwr0xdTFTy07NdxeR5Cj4cK+tKueK2kpWK9xFJEEKPtwh1Xtv2XWUnh7dRExEkkHhTmrc/dipTt5rOxl1KSIioVC4A/MbxwJoSqSIJIbCHWgcV0FtVZnG3UUkMRTugJnxiSvG8eZ7h7V4h4gkgsI9cNOMcXxw/Aw7DnVEXYqIyJAp3AM3zagFYOV7hyOuRERk6BTugcZxFUyqHsGbCncRSQCFe8DMuHHGON58/7Dmu4tI7Cnce7l5Ri1HOs6x5YMTUZciIjIkCvdebpwxDoCV78VirRERkYwU7r1MrhnJ9NpKjbuLSOwp3Pu4ccY4Vu04Qld3T9SliIgMmsK9j5tn1HLybBcb9rVHXYqIyKDlskD2D83soJlt7KfNIjNbZ2abzOxX4ZZ4eS24InWfGc13F5E4y6XnvhS4K9OLZlYDPAV8zt2vAT4fTmnRGFdVzlUTR/HGuzqpKiLxlTXc3f11oL87an0BWObuu4P2B0OqLTILm+pYvfMIJ850Rl2KiMighDHmPhMYY2avmdkaM3swU0MzW2xmLWbW0tbWFsKm8+P2qybQ1eP8i3rvIhJTYYR7CXA98GngTuCPzGxmuobuvsTdm929ua6uLoRN58e8hhqqR5ayvDX2/wgRkQJVEsJn7AUOuXsH0GFmrwOzgW0hfHYkSoqLWNRUx2tbD9Ld4xQXWdQliYgMSBg995eAW8ysxMwqgE8ArSF8bqRuu2o8hzvO8c7eY1GXIiIyYFl77mb2PLAIqDWzvcA3gFIAd3/a3VvN7OfAeqAHeMbdM06bjIuFM+soLjJebT3IvIYxUZcjIjIgWcPd3e/Poc23gW+HUtEwUVNRxvXTxrB8y0H+y51NUZcjIjIgukK1H7dfNZ7WA8fZf+x01KWIiAyIwr0ft189HoBXt2jWjIjEi8K9HzPqqmgYW8Hy1g+jLkVEZEAU7v0wM+6YNYH/t/0w7ad0taqIxIfCPYvPzZnMue4efr7pQNSliIjkTOGexbVTqpleW8lL6/ZHXYqISM4U7lmYGZ+dPZk33z/MweNnoi5HRCQnCvccfG72ZNzhn9ZraEZE4kHhnoOPja/imsmjeXndvqhLERHJicI9R/fOmcw7e9vZcagj6lJERLJSuOfos7MnYwYv68SqiMSAwj1Hk6pHMr9xLC+t24e7R12OiEi/FO4D8FvN9bx/qEOLZ4vIsKdwH4DPXDeJsZVlLF25M+pSRET6pXAfgBGlxfzO/HqWt37IniOnoi5HRCQjhfsA/bsF0wD4u1W7Iq5ERCQzhfsATa4ZyR2zJvL3q/dwprM76nJERNJSuA/Cl25q5NipTk2LFJFhK2u4m9kPzeygmfW7LqqZzTezbjO7L7zyhqcFV4ylacIolq7cqWmRIjIs5dJzXwrc1V8DMysGvgX8IoSahj0z46FbprP5wHF+uVkLeYjI8JM13N39deBIlma/B7wIFMx6dL85dwpX1FbynX/eSnePeu8iMrwMeczdzKYAvwE8PfRy4qOkuIj/fMdMtn14kp+8rRuKicjwEsYJ1e8BX3X3rFNHzGyxmbWYWUtbW1sIm47WPR+fxDWTR/Pd/7uNc109UZcjInJBGOHeDPzYzHYC9wFPmdmvp2vo7kvcvdndm+vq6kLYdLSKioz/emcTe4+e5vl/3R11OSIiFww53N19urs3unsj8ALwu+7+kyFXFhMLZ9Zxw/Sx/M9Xt9N+Wotoi8jwkMtUyOeBN4EmM9trZg+Z2SNm9kj+yxv+zIz/8emrOdJxlj95ZXPU5YiIAFCSrYG735/rh7n7l4dUTUxdN7WGRxbO4KnX3uPuaydy21UToi5JRAqcrlANyWOfupKZE6r4w2UbaD+l4RkRiZbCPSTlJcX82efncOjkOf74lU1RlyMiBU7hHqJrp1bz6KIZLFu7j797S3eNFJHoKNxD9vu3X8mtTXV8/aWNLG/VrQlEJBoK95CVFBfx/S/M45rJ1fzH//U26/cei7okESlACvc8qCwv4dkvNzOuqox/v3Q1rQeOR12SiBQYhXuejB81gqVfuYHiIuO+H6zk1S0aohGRy0fhnkcfG1/FS49+kul1lTz8oxaefWOH7v8uIpeFwj3PJlaP4H//hxu5Y9ZEvvnKZr7016vZcagj6rJEJOEU7pdBRVkJT31xHt/47CzW7jrKnd99nSf/eSsnzuhiJxHJD4tqmKC5udlbWloi2XaUDh4/w5/+tJWX1u2nsqyYf3v9VB68cRofGz8q6tJEJAbMbI27N2dtp3CPxjt7jvGjN3fyyjsHONfdw9WTRrNwZh0LZ9Yxp76GkWXFUZcoIsOQwj0mDp08y7K1e3l1y0Fadh6lq8cxg+njKmmaOIpp4yqZVD2CidUjqK0qY/SIUkaNKKWivJgRJcWUFhtmFvVuiMhlonCPoZNnu3jrvcNs3N/OlgMn2PLBcfYdO01nd+ZjZAalxUWUFBnFwY+RuhWxBa+DXWib9jP6fF7qOev1Ox/5ArnwvKXaXfpeLvnCsQwPPrrtwX9J6etN4uS359fz8C1XDOq9uYZ71lv+yuVTVV7Cp2ZN4FOzLt4yuKfHOdxxjgPtpzl6qpPjpzs5fqaT0+e6OdvVw9nObs5299DT43T1ON09jjs45/9MyfwdfvGF823Ov//i75c+j/d+JxemeHqabaVr1/d5htDH8KG8WSQCtVXled+Gwn2YKyoy6kaVUzcq/38ZRCQ5NBVSRCSBFO4iIgmkcBcRSaBcFsj+oZkdNLONGV7/opmtD35Wmtns8MsUEZGByKXnvhS4q5/XdwAL3f064JvAkhDqEhGRIcg6W8bdXzezxn5eX9nr4VvA1KGXJSIiQxH2mPtDwM9C/kwRERmg0Oa5m9mtpML9k/20WQwsBmhoaAhr0yIi0kdOtx8IhmVecfePZ3j9OuAfgbvdfVtOGzZrA3blXOlH1QKHBvneOCvE/S7EfYbC3O9C3GcY+H5Pc/e6bI2G3HM3swZgGfBArsEOkEtx/WyzJZd7KyRNIe53Ie4zFOZ+F+I+Q/72O2u4m9nzwCKg1sz2At8ASgHc/Wng68A44Kngxk9dhXiARESGk1xmy9yf5fWHgYdDq0hERIYsrleoFupc+kLc70LcZyjM/S7EfYY87Xdk93MXEZH8iWvPXURE+hG7cDezu8xsq5ltN7M/iLqefDCzejNbYWatZrbJzB4Lnh9rZr80s3eDP8dEXWs+mFmxmb1tZq8Ej6eb2apgv//ezMqirjFMZlZjZi+Y2ZbgmN9YCMfazP5T8Pd7o5k9b2Yjknis092fK9PxtZS/CPJtvZnNG+x2YxXuZlYM/CVwNzALuN/MZkVbVV50AU+4+9XAAuDRYD//AFju7lcCy4PHSfQY0Nrr8beA7wb7fZTUxXJJ8ufAz939KmA2qX1P9LE2synA7wPNwfUzxcDvkMxjvZRL78+V6fjeDVwZ/CwGfjDYjcYq3IEbgO3u/r67nwN+DNwbcU2hc/cD7r42+P0Eqf/Zp5Da1x8FzX4E/Ho0FeaPmU0FPg08Ezw24DbghaBJovbbzEYDvwY8C+Du59z9GAVwrEnN1htpZiVABXCABB5rd38dONLn6UzH917gbzzlLaDGzCYNZrtxC/cpwJ5ej/cGzyVWcHXwXGAVMMHdD0DqCwAYH11lefM94L8BPcHjccAxd+8KHiftmF8BtAF/HQxFPWNmlST8WLv7PuA7wG5Sod4OrCHZx7q3TMc3tIyLW7inW+Q+sdN9zKwKeBF43N2PR11PvpnZZ4CD7r6m99NpmibpmJcA84AfuPtcoIOEDcGkE4wx3wtMByYDlaSGJPpK0rHORWh/3+MW7nuB+l6PpwL7I6olr8yslFSwP+fuy4KnPzz/T7Tgz4NR1ZcnNwOfM7OdpIbcbiPVk68J/ukOyTvme4G97r4qePwCqbBP+rH+FLDD3dvcvZPULUxuItnHurdMxze0jItbuK8GrgzOqJeROgHzcsQ1hS4YZ34WaHX3J3u99DLwpeD3LwEvXe7a8snd/9Ddp7p7I6lj+6q7fxFYAdwXNEvUfrv7B8AeM2sKnrod2EzCjzWp4ZgFZlYR/H0/v9+JPdZ9ZDq+LwMPBrNmFgDt54dvBszdY/UD3ANsA94D/nvU9eRpHz9J6p9i64F1wc89pMaflwPvBn+OjbrWPP43WETqTqSQGpf+V2A78A9AedT1hbyvc4CW4Hj/BBhTCMca+GNgC7AR+FugPInHGnie1HmFTlI984cyHV9SwzJ/GeTbBlKziQa1XV2hKiKSQHEblhERkRwo3EVEEkjhLiKSQAp3EZEEUriLiCSQwl1EJIEU7iIiCaRwFxFJoP8PYxnFT1B9cX4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10894e080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXd4Y9d55t+DXggQIACSw16mz2iaR2VULFmSHcmyLa9TbMeJ7dhZrb1O4mzyxCXr7Mab3eymrOPESZwo7o5jO+siy5KbeplRG03vQ3IKOwGQAIjezv5x7wVBEJW4AIiL7/c880gkwYsDXuC9333PVxjnHARBEIRyUDV6AQRBEIS8kLATBEEoDBJ2giAIhUHCThAEoTBI2AmCIBQGCTtBEITCIGEnCIJQGCTsBEEQCoOEnSAIQmFoGvGkTqeTDw0NNeKpCYIgmpbXX3/dwzl3lXpcQ4R9aGgIR48ebcRTEwRBNC2MsWvlPI6sGIIgCIVBwk4QBKEwSNgJgiAUBgk7QRCEwiBhJwiCUBgk7ARBEAqDhJ0gCEJhtIywp9Mc//7aJOLJdKOXQhAEUVNaRthPTvnwie+fwvOX3I1eCkEQRE1pGWFfjiYBAL5IosErIQiCqC0tI+zheAoAECBhJwhC4bSQsAsReyBKwk4QhLJpGWEPZSL2ZINXQhAEUVtaRtgjYsTuJyuGIAiF0zLCHoqJETtZMQRBKJyWEfaMx04RO0EQCqdqYWeMGRhjrzLGTjLGzjLGPivHwuQmkxUTJY+dIAhlI8cEpRiAuznnQcaYFsCLjLGfcs5fluHYskHpjgRBtApVCzvnnAMIil9qxX+82uPKTShGVgxBEK2BLB47Y0zNGDsBYAHAE5zzV/I85iHG2FHG2FG3u/5l/ZGEELEvx5JIpTfcdYcgCEI2ZBF2znmKc74PQB+Amxhju/M85mHO+UHO+UGXq+SQbdmRInYACJLPThCEgpE1K4Zz7gPwLID75DyuHEgeO0ApjwRBKBs5smJcjDGb+P9GAPcCuFDtceUmHE/BpFMDoCIlgiCUjRwR+yYAzzDGTgF4DYLH/pgMx5WVcDyJ7nYDANpAJQhC2ciRFXMKwH4Z1lJTwvEUtnVbMOEOkRVDEISiaYnK03SaIxxPodtqBECNwAiCUDYtIexSquMmyYqhiJ0gCAXTEsIuZcS4LHqoGG2eEgShbFpE2AXrpU2vgcWgpc1TgiAUTUsIu9Sy16RTw2rUUCMwgiAUTUsIeyQhCLlJr4GVInaCIBROSwh7dsTebtTS5ilBEIqmJYRd8thNOjWsBi1tnhIEoWhaRNiFiN2s0wgeO+WxEwShYFpC2EPxrM1TA1kxBEEom5YQ9nBsZfO03ahFOJ5CIpVu8KoIgiBqQ2sIuxixG7VqWI1aANQIjCAI5dIiwp6EQauCWsVgNQp9zyiXnSAIpdIiwp6CWScIutVAETtBEMqmZYTdpBeGbGSsGNpAJQhCobSEsIdiSZi0QsTenvHYyYohCEKZtISwRxJZEbtoxVCREkEQSqUlhD0US6547JnNUxJ2giCUSUsIezieglEcZG3UqqFRMdo8JQhCsbSMsJtFYWeMwUqNwAiiZRl3BxGKKXuPrUWEPQmTfmVud7tRS5unhGx8/chV/PIXj4Bz3uilECWIJlJ4+xdexGceOdPopdSUlhD2UCwFk1ad+dpq0NDmKSEbPz0zi9evLeHsTKDRS2l5Ls0vF73Avn5tCeF4Co+enMHkYjjzfc45/st3T+Afnx2rxzJrjuKFPZ3mYlbMSsROVgwhF6k0x6kpPwDgqfMLDV5NazPuDuItf/M8nrlY+DwcHvNArWJQMeBLL0xkvv/TM3P44fFpPHvBXY+l1hzFC3skIbXszY7YaYoSIQ+X5pcRjqegVjE8dWG+0ctpaWZ9UQDA+dnlgo85PO7Fvn4bHtzXi+8enYQ3GEM4nsT/fOwcAGBhOVqXtdYaxQt7KGvIhoQQsZPHTlTP8es+AMCvHOjDqSk/5gPKEIZmZCkcByBE7vkIRBM4PeXDbaMOfOTOEUQTaXz9yFV88dlxzPij2D9gw3wgpoi9EsULeyTTiz3biiGPnZCHE5NL6DDr8Fu3DwEAnr5Adkyj8InCfsUTyvvzVyYWkebArZud2NxpwVt2duGrR67in5+fwDv39eD+3d2IJFJYVkDGjOKFXZp3atavtmLiyTSiok1DEOvl+HUf9vfbsK3Lgl6bEU+eIzumUSyFhWBtwh3KG3UfHvPAoFVh/4ANAPCRu0axHE1Cq2L49Ft3oMtqAAAsKOCuq2phZ4z1M8aeYYydZ4ydZYx9XI6FyYU079SoW715ClD1KVEdgWgCY+4g9vXbwBjDvTs68eKYJ3OXSNQXyYrxRxIZkc/mpXEvbhzqgF4jBHkHBuz48O3D+LN37kaX1YBOiyTssfotukbIEbEnAfwh53wHgFsAfIwxtlOG48rCyrzT1emOADUCI6rj1KQfnAP7B+wAgHt2dCGWTOPwmGfNYyPxFH56elYR/u1GxZcl5hM5Prt7OYaL88u4ddS56vt/8radeNeBPgBAp1UPAJhXwAZq1cLOOZ/lnB8T/38ZwHkAvdUeVy7Cmc3T1QVKAEXsRHUcv74ExoA9/e0AgJtHOmDWqfNmx/zNk5fw0W8do1z3GuILxzNB20SOz35kXLjY3rbZUfD3JStmniL21TDGhgDsB/BKnp89xBg7yhg76nbXL1c0nDXIWkKyYmgDlaiG45M+bHa1ZTqG6jVqvHGrC0+eX1i1fzPnj+LrR64CAM7PkrDXiqVwArt726FVM0y4Vwv7S+NeWA0a7OppL/j7bXoNzDo1WTHZMMbaAHwfwO9zzte8eznnD3POD3LOD7pcLrmetiQhSdhzNk8BmqJErB/OOU5M+jIbcRLvu3kQ7uUY/vwn5zPf+8LTl5FKc+jUKlyaL5xjTVSHLxyHo02PQYd5jRVzeNyDW0YcUKtY0WN0WQ1kxUgwxrQQRP1bnPMfyHFMuQiLqUvmnHRHgOaeEuvn+mIYi6E49vXbV33/9i1OfPj2YXzjpWv4+dk5XPOG8N3XJvHemwawpasNF+ZI2GvFUjgBu0mLYad5VcrjVU8Ik4sR3LbZWeS3BVwWPWXFAABjjAH4MoDznPPPVb8keZGsGKOWInZCPqTCpNyIHQA+cd827O614hPfO4X/9qOz0KgZfvfuzdjWbaGIvUak0hyBaAI2kw4jLjOuecNIpYWNaqnFwJu2dZY8TpfVgIVlsmIA4DYAvwngbsbYCfHfW2U4riyE40kYtWqosm7BDFo19BoVCTuxbk5M+mDSqbG1y7LmZ3qNGl947wEkU2k8d8mND946jE6rAdu6LJgPxDKFNIR8+CMJcA7YTVqMOM2Ip9KYXooAAJ656MaIy4wBh6nkcbqseswHok2fvSRHVsyLnHPGOd/DOd8n/vuJHIuTg1A8tWrjVKLdqKXNU2LdTC6GMeQwF/Rsh51m/PWv7sWNQ3Z85M4RAMC2buEicJHsGNmRLpY2kxYjrjYAwLgniHA8iZcnvGVF6wDQaTEgmkg3vU2rKf2Q5iYST63aOJWgDo9ENXiCMbgs+qKPuf+GTbj/hk2ZryVhvzS/jJtHCqfdEZUjFSTZTDoMO80AgCvuENJpjngyXb6wi7nsC4FoJi26GWmBlgLJVRunEjRsg6gG93IMzrbiwp5Lt9UAi0FDG6g1QIrY7SYdHGYdrAYNJjxBPHNxASadGjcO20scQSDTVqDJfXbFR+zZ806zsRo08ATJ6yQqh3MOTzAOp0VX0e8xxrCdNlBrghSx201aMMYw7GrDhDuEa94wbtvszLQRKMVKkVJzZ8YoPmIPx/NH7GTFEOslEE0inkrDVWHEDgBbuyy4OFd8yg9ROSseu3CxHXWacfTaEqZ9kbJtGADoFO21Zq8+bQFhp81TQl48QeFDX6kVAwg+eyCaxFyTR4QbDV84AbWKZVoKjLjMiCfTAIC7tpVfEGnWa9Cm1zT9wA3FC3sonswr7NIUJYqciErxLFch7F2UGVMLlsJxtBsFGwYAhp1CZsz2bgt6bMaKjtVp1Td9WwHFC7uQFZN/8zTNgaACmuoT9cUtRuylsmLysZWEvSb4wgnYTCtZLCMuITPmrgpsGIlOi5489o2GsLG1crUNxVKrWvZKUFsBYr2sROyVbZ4CgN2sQ6dFj4u0gSorS+E47KaV87G1y4Lfu3sz3n9osOJjKaFfjOKE/ZmLC7jlz5/C6Sk/0mmOSCK1asiGhJSj6s/TkJ8giuEJxqFWsVVCUgnUWkB+pD4xEmoVwx+8ZVvFNgwgthVo8tmnihP2aV8UyTTH3z19GZHE2iEbEpl+MZQZQ1SIJxhDh1m3qk1FJWzrsuDyfDDTy4SoHl84nsmIqZZOix6xZLqp61wUJ+xSN8cnzs3j6LUlAMjrsWfG41FmDFEhnmDlxUnZbO22IJZM46o3/9BlonJ84QRsMlWKdkq57E1sxyhO2EOisFv0GvzVzy8AAEza/OmOAA3bICrHvVy6nUAx9vcLHSFfnvDKtaSWJppIIZJIwW6WJ2LvskhtBZo3M0Z5wh4XNkt/67YhnJkW5n2Y8/WKyVgxzXu7RTQGTzC+ro1Tic2dbRh0mPDkubUj9IjK8WX6xMgTsSuh+lR5wh5LwqTX4EO3D2e8dVOezVOLQQPGlBmx//XPL+Kfnhtv9DIUCecc7mBsXVWnEowx3LujC4fHvZk7TGL9LGX1iZEDJQy1Vp6wx1No02tgM+nw/luHAOSP2FUqhja9RpEe+2OnZvDPz40jmUo3eimKYzmWRDyZrspjB4B7d3QhnkzjhcsemVbWuixlteyVA5NOA4teQ1bMRiIUW6k0/c93jeKP37ode/rWTrkBpA6PyhN2bzCOpXAis3lMyEcmh73CBmC5HByyo92oxZPnyY6pFl+mAZg8ETsgVp9SxL5xCMWSMItZMBaDFg+9cRRadf6XaTUorxFYNJHCsnh7/4uzJBpyI3UEdbUZqjqOVq3Cm7a58PSFhUzaYzSRwt88cQlz/uYVlEYgt8cOAN3tBkz7mvc8KE/Y48m8eev5sBo1TZ2rmg9vSBAeFQN+cW6uqYssNiJumSJ2ALh3ZxcWQ3Ecvy7cWX32x2fxt09dxmOnZqo+disht8cOACPONky4g037+VGcsIdjqUzEXgoldnj0iu0U7t7ehamlCM7PUoWjnFTT2TGXN251QatmeOL8PL73+hS+/eokAOCaN1z1sVsJXzgOg1YFQ5605vUy6jJjOZrMXMibDcUN2ggWmJiUDyVaMV7RKnjPjf146sI8fnFuDjt7rA1elXLwBGNQMXmiQ6tBi1tGHHjk+DT8kQRuGenAcjRJhUsVIrQTkC9aB4DNnUKztjF3MFOw1EwoL2KPt3bELkWU27otODhoJ59dZoR2AvqCQ6wr5d4dXZgPxGA1aPGF9x7AqKuNhL1C5GwnIDHaKXSHHF8IynrceqEoYeecCx57nvTGfFiNWoTjKSQUlBYoeeyONh3esrMb52YDmFykW3u5cC9XV5yUy/03dGNfvw3/+L4DcFn0GHKYML0UyQyJIEqzJGM7AYluqwFmnRrj7ua8yCpK2COJFDhHRRE7oKx+MZ7lGIxaNUw6Dd6yqwsA8AuqcJQNd7C6dgK5dFoMeORjt+HgUAcAYMhpRpoDk0t0MS4XXzgOu1leYWeMYbSzDWMUsTceaWhGJVkxgLLaCnhDcTjEiHLQYca2LguevkDCLhee5eqqTksx6BAsgGtkx5SNMGRDXisGADa72jDuJmFvOOGY2Ka3wohdST57bufB3b3tGF8gkZADaYiLU8aIPZchhwkAcMVDEXs5cM7hi6zuxS4Xo51tmPVHy5qy9uJlD971j4cRS6ZkX8d6UJSwSycgX2+YfGQagSlI2L05DaoGHSbMBaKIJjbGG66ZCcaSiCXTsnrsuXSYdbAYNBSxl0kgmkQqzWXPigGAUZcwN7WcDdRnLy7g2HXfhhl5qChhD8cF8WorM2LP9GRXUMqjNxSDw7wSUQ50CBHgFHm2VSNVncqRw14IxhiGHGZc8ZCwl4NPLE5ql3nzFAA2S5kxZdgxUibT2ZmA7OtYD7IIO2PsK4yxBcbYGTmOt16kTnmmMrNilGbFpNMc3uCKxw4AA+KtPRW9VI9UrCLn5mk+Bh0mOl8AUmle0gZZqkGfGIlBhxkaFStrA1W6EJ9TkrAD+BqA+2Q61roJxYU3QdkRe8aKUcbmaSCaQDLN4WhbG7Ffp5THqpGz6rQYw04zppbCLZ/y+O1Xr+P2v3i6qI14ctIHAOjrqHy2aSm0ahUGHaaSEXsqzTG5GAEAnJ3xy76O9SCLsHPOnwewKMexqiETsZeZFWPQqqBTqxQTsa9YBSvRi8Osg1mnpghQBuol7IMOIeWx1e2zsYUgfOFEQd+ac45vvnwNe/rasa3LUpM1jLpKpzzO+CKIp9KwGjQ4P7u8IWbZ1s1jZ4w9xBg7yhg76na7a/IcoVhlHjtjTGgEphCP3ZtHeBhj6O8wUZGSDHiWhXYCHTKNYCvEsJPsMwBYFIvtCvnWL417MbYQxPsPDYExeSqBc9nc2YZr3nDRIkbJhnnLrm5EEqkNUTlcN2HnnD/MOT/IOT/ocrlq8hyhCrNiAGEDVSkRe3bVaTaDDhOukbBXjTsYR4dZJ1s7gUJIueytvoHqDQmBSiF74+svXYXdpMXb9myq2RpGXW1IpnnRi6wk5A+I69gIG6iKyooJxVPQqVXQacp/WVaDcoZtSFZBdlYMIAjF5GIY6Q1wi9jMuJdjNbdhAME+a9NTyqPU0C6fUE77Inji3DzefeOArF0dc9ncKaQ8FrNjrnhCMOnUuHXUAa2abQifXVnCHkuWnREjYTVqFVN56gnGwRjWFGv0d5gQS6ax0KQtSDcKF+cDmdzmWsIYw5DThKtlWDEnJn0YW9gYudNyI92Bnp8NrBnz+G+vXAMAvO/mgZquYcRVOuXxqieEQYcZeo0aWzotGyIzRq50x28DeAnANsbYFGPsw3Ict1KEIRuVdSJW0ng8bzAGu0kHTc7EqMEOybNt7QiwGpZCcUwuRnBDX3tdnm/QYS7Lq/34d47jI/96THF3Y+k0x2Iojp52A2LJNCaybKlYMoXvvDqJe3Z0oV98b9cKi0GLbquhaJHSNW84sy+yq8eKczOBhg/okCsr5r2c802ccy3nvI9z/mU5jlspwli8CiN2g3IGWudWnUpQymP1nJ4Wbq/39NZH2IcdZkwtRYpu2kUTKVxfDGNsIYjnLtUmIaFRBKIJpNIcd2wR9uOy7Y2fn52HNxTH+w8N1mUto51mXCiQmZNMpXF9MYwhcV9kV48V3lAc8+Ig7HF3EP/xG0frnqChKCumkl7sElJP9kZfYeXAE4yt8dcBoNduhIqRsFeDJOy76iTsgw4TUmmOqaVIwcdMuEOQ3rb/8sJEXdZVL6TU3ZuGO6DXqHB2esXe+PHJGXRZ9bht1FmXtdy1tRPnZgM4M73WO5/2RZBMcww5BWHf2SO8P87N+sE5x6e/fxpPnJvH5fn62mWKEvZKpidJWI1aJNMcEQX0Usnu7JiNVq1Cj81Iwl4Fp6Z8GHaaa1K6ng9JKIrZMZLv+859PTgy7t0Qm3ZyIaU6dlr12N5tyWyg+iMJPHfRjbft6YGqxtlJEr92Yz9MOjW+cvjKmp9JmUtSxL5jk5BPf3Y6gO8fm8arV4Xynki8vsVmihJ2Yd5pZVaMktoK5HZ2zIbK1Kvj1JQfN9QpWgdWhOJKkUEPYwtBqBjwxw/sgEmnxpdfWCs8zYo3K8NrV287zs4IEfAT5+YRT6VrmuKYS7tRi1872I8fn5zBQiC66mdXJWEXPXaLQYshhwlHxr343z85n/k81jtwVJSwrytiV0hbgVgyheVoEo4CxTMDVKS0bhaWo5j1R7GnThungFA97GzT40yRKHzcHUR/hwmdFgN+7WA/Hj05gzl/tODjm4nsmoxdPVYEoklMLUXw45Mz6LMbsa/fVtf1fPDWISTTQqVrNle9YZh16lU9+nf2WPHShBe+SAL//e07AZCwV0U4nlyXxw40f8Qu3boW6hU+0GGGNxQvq7c0sRrJW93TVz8xYYxhb187Tk0VFvaxhWAm/fLDtw8jzTm+duRqnVZYW6QcdrtJh12ib/3imAeHxzx4256emlWaFmLIacY927vwrVeur+pdc8UTwpDTvGo90no/eOsQ9g8I75l6t81WlLCHYql15LGLU5SaXNg9y2KEUyRiB4DrZMdUzKkpPxgTMh7qyd5+G8bdQSznyahIpTkmPKFMAU1/hwlv2taJx0/P1HWNtWIxFIPVoIFOo8L2bgvUKoa/f3oMyTTH2/fWz4bJ5kO3D2ExFMcjx6cz37vqDWX2QyQeuGETfu1gH37/3i0wisVTJOzrJJ5MI55Ko229VkyT94vxiOXXjiIeOwBcX1zr2R4Z9+CT3zuliMygWnB6yo/NrraK7warZU9fOzhfycjJRhp4PepaEZU7tjgxuRhRhOXmCcUz/rRBq8aoy4xpXwQjTjN2bqrvBVbi0IgD27st+KfnxuGPJJBIpTG1FMGwY7WwDznN+Mtf2QuLQQuj2JAwEidhXxfhuNSLvTWtGG+ezo7Z9BfJZf/Ki1fw3aOTGTuHWIFzjpNT/roVJmWzV7R+Tk6uFfYxt5A+J0XsAHDbZiH97/CYpw6rqy3eYGxVhpdkb7xtb/1tGAnGGP7kbTsx7Yvgt7/+Gi7NC50ccyP2bAwaUdgpYl8focz0pMpbCgBKEPbiEXu7UQubSbsmMyaaSOFFUQg2ajrki5c9+NNHzzbkuecCUXiCsboVJmVjN+sw0GHCqSnfmp9Jc2yzWxxs7mxDp0WPw+Peuq2xViyG4qu6aO4VL6xvr2M2TD5u2+zE59+9H69fW8JD33gdwEo3znyoVAx6jYqEfb2sp7MjAKhVDBa9Br5wkwt7KA69RgVzkV70Ax2mNeL90rgX0YSQY7tRhf1fX76Grx252pCNX2nzck+dszAk9vbbMsMkshlbCMLZpoMta3IQYwy3jjrw0rin6W01YRLYSpDynpsG8KOP3YYtNeq7XgkP7NmEP/8PN2DaJxSPDTkKR+yAYCXFEpTHvi4kYS+3F3s2VgX0i/GInQeL3aaOOM04OxNYNUn96QsLMGiFt8FG3FjlnOPY9SUAjel1c3rKD7WKNczX3dvXjhl/NDOWT2LcHcRInoZkt252whOM42KdKx3lJJXmWArHVyUCGLRq7G3QxTUf77lpAJ99xy7cvb2zZH9+o1ZNHvt6kYZslDs9KRubqfl7sgubTcXfYO860IfFUByPnhAyJzjnePrCAu7Y4kKXVb8he7ZP+yKZrpSNKLA6Ne3H1i5LTVvDFkMSs2w7hnOOMXdwlb8useKzN68d4wvHkeaFM7w2Ch+4dQhf+eCNJT1/o05NVsx6keadridzoV0BwzaEzabivcLv2OLE9m4LvvziFXDOcXF+GdO+CO7Z3pnXptkIHLu+ImiNmEwzvRQu6qHWml09VqgYVtkxi6E4fOEENueJ2HttRqHysYk3UKVN/I469L6vBwYtCfu6kayY9Qq7r8mFfT4QXVX9lg/GGD58+zAuzC3jhcsePHV+AQDwpu2dG3Z83rFrSzBoVegw6xpiFXlzNvHqjUmnwdYuC05mFSpJQx9G80TsgGDHvHJlcU0P82YhM7t3g0fs5WLQqiiPfb1IWTGV9ooBmt+KCUQT8ATjGHYV38QBgHfs64HLose/vDCBZy4sYHevFV1WAwY7zJgLROv+BizF8etL2NNrw4izvP7kcpJMpeELJ9CRp2NmPdnbZ8PJKV9mQ3RMbP6Vz4oBgNtGnQjGkqsuBs2Et0RNRrNh1KpJ2NdLJmKvMCsGWJl72qyZBFKjqOEi+bQSeo0aHzg0iBcue/D69SXcvb0LADDgMIJzFG0TW2+iiRTOzgSwf9CGgQY0MVsSM6VK7V3Umr39NvjCCUwuCudmfCEEo1aNTVZD3scfGnUAQNPaMRkrRiERu5GsmPUTjiXBGDIlvJXQbtQinkxn0v6aDal16EgZwg4A77t5EAatCpwD92zvBCD0kgGwoeyYM9N+JNMcBwbsGHKYMeuv7x3FRhEYqfnYY6dnMLYQxPnZAEY7zQXb1naYddi5yYrD480p7IVGPDYrBh1lxaybYCwFk1a9rh7NNqPwwW1WO2bCEwJjwICjvE0+u1mH9x8awpDDlGlFO7ABx+dJaY4HBuyZlgj1vPBIlkCjhX1btwUWgwZ/+bOLuPdzz+GlCS+2dBbP596xydq0bZoXQzHYjNo1Ix6bFcGKqW/QWN/mFzVkPZ0dJbLbCnS357+93chc8YTQZzdCryn/buVT923HH/3StsyF0Nmmg0mnxvXFjWPFHLvmQ3+HES6LPlMEctUbrluRihSx55tKVU+0ahV++vE7MO4OwReOIxBJ4M6tnUV/p9Oqh3s5hnSa120ghVzkFic1O43YPFWMsAdj1Qu7L9ycvVKueIIYdubfSCuESsWgwsoHnjEmpjxujIhdKky6VfSLpYg9+47iujeMrx65gk/et70meeYbxYoBgD67CX328tMuOy16JKVCnyYTyUZnIskNeexVIMw7Xd+H22Zq3n4xnHNccYfK9teL0b+BctmlwqQDg3YAgM2kQ7tRuyoz5luvXMNXD1/FN1+6VugwVbHSE7z5vN5Oi3Dn6Q7GSjxy4+ENxhq+YS0nkrDXMzlDMcIejCUr7hMj0cwdHt3LMYTiqbIyYkohFSlthOwgqTDpwIA9872hnMyY5y65AQB//8xYTc6dNxSDzdScXq9LHLiyEGg+YV8MxRtuf8mJQacG50AsWT+fvfnesQUIx5Pr6hMD1KbD43wgWpd0s3Ex1XGkjBz2Ugw6TIgm0mv6kjSCE9d9MGiFIQsSgw5zRtjn/FFcmFvGO/b2wB9J4J+fG5d9DYLANGfk2CkJ+wY4l5WQTKWxFE4ozooB6jtsQznCHkutq08MAFj0GqiYvML+xWfH8RtffqXm1ZJSqqMcEXuxnu31Zi4QQZ/dtCpaHnSYMO2LIJFK43kxWv/oXaN4cF8PvnL4CuYD8s779AabN3LWKgDBAAAgAElEQVTMROzLzTUDdTFcfK5AM2LICDtF7BUTjK0/YlepWKZISS4uLywjzYGvHqnt5PgrniB0GhV62o1VH2twAwn7YiiODtPqD/egw4xUmmN6KYLnLrnRZdVje7cFf/jmbUilOT7/5GX519CkkaNZr4FZp94Qd1+VsLJh3ZwX1HxIEXs9N1AVI+zheGrdHjsg9ouRsSe7NAjh31+brKl3f8UTwrCjcLFKJfTajWCsMV0Uc1kMxWE3r960HBIzY8bdQbxw2Y07t7qEbB6HCe+7eRD/fnQSMz750jUXQ3F0NHHk2Gk1NJ0VI21YO5r4756LFLHXs0hJEcLOOUconqx4elI2Nhkj9lAsiblAFG+9oRuheArfefV6yd+ZD0Txvi+9jIUK7YQJT0gWGwYQ2g1sshoqKgLinONrh6/AL/OgksXQWp91UMxlf+TEDALR5Kpc7ncd6EUqzfNOG1oP6Tw9wZsNl0UPd5NtnnpDxYeyNyOZuafNFrEzxu5jjF1kjI0xxj4lxzHz8cVnx/GbX35lzfeFVKLK551mI6cVI/neb9/Tg0MjDnztyFUkSnTae+r8Ag6PeTPVluWQTKVx3Rsuq/lXuQw4TBX1Zb/iCeFPf3wOj56akW0NnAuias+xYqQiqp+enoWKAbeLvccBoSEWY8DFuaAsa/BFEkjzjZHDvl46LfqmS3csNeKxGWnKzVPGmBrAPwC4H8BOAO9ljO2s9rj5CMWSODLuXSOSwSpa9kq0yzhFady90lb1t+8Yxqw/ip+cni36O5Kgz1cQYU0tRZBMc9kidiD/+LxiSO2O5bRAAtEkUmm+RlQZYxh0mJFMc+wfsKM9K7/cpNNgoMOESwvyTA5aVECHQZdFX/EdYKNZDMWhYsIdtFKQJpQ1lbADuAnAGOd8gnMeB/AdAA/KcNw1DDlXNs+yCYvTk4rN+yyFnD3Zx90hqJiQxfGmbZ0YcZnxLy9MFM0Pl4R9roIPYqXNv8phoMME93KsbD/QXwNhXypS8Sn57Hduda352dYuCy7NySPsUk/wZrYEOi0GhOKpTOfTZuDklB9DTnn2jDYKzbp52gtgMuvrKfF7q2CMPcQYO8oYO+p2u9f1RNKH+kpOoyo5InapJ7scxTkT7iD67CboNUJTsncf7MeZ6UBGLHLxheOYEPPRK0nZm5Ax1VGiW8yuKTdNTrrLmfXJFxlKKW/2PKIq+ez5hb0NVzyhVTNd172GDdROYL1IuezNkhkTjCXx8rg303FUKTTr5mm+S+sadeScP8w5P8g5P+hyrf1QloP0ob7mWS3sYWnIRpVZMak0zwzsqIZxdwijWb73VrFpVaFBEcfFKkudWlVRpeCEOwirQSOr+Ej5w54yvVlJ2KdrEbGb1r6ud+ztwW/dNpTpSpnN1i4LkmmeuZOpBiVs4rmarEjpxcsexFNp3LOjq9FLkRVp87TZrJgpAP1ZX/cBkG8nLQtnmw5mnRpXc9LxVsbiVWfFANU3AkunOa54Vk+QHxIj6qsFBOfY9SWoVQw3j3RUFLFf8YQw4morOUy3EpxtUpRX3t9BsmLmA1Gk0vK0IvAWiZZ39ljx39++K++t+jaxSvXSfPUbqIvBwncNzUKntbmKlJ46Pw+rQYM3DNpLP7iJMDZpgdJrALYwxoYZYzoA7wHwqAzHXQNjDEN5RqRVM8haol2mnuwz/giiifSqEv8+uxFqFSsYsR+7voTt3RaMOM2VC7uMNgywIuzlRuzS3yuZ5mX/TimkiL1SUR1xtkGjYrL47IuhGKwGDbRN2CdGQmoE1gz9YtJpjmcuLuDObZ1N/TfPh6EZPXbOeRLA7wD4OYDzAP6dc3622uMWYiirX4hENYOsJeRqBCZ55aNZEbtWrUK/3YirnrXZJqk0x4nrPhwYsKPTakAgmizLi1uOJjDrj8rSIyYbR4VWTPbfSy47ZjEch06tqngzXKdRYchpxsX56oXdG2q+dre52IxaaNWsKVIeT0754AnGce8OZfnrAKBWMejUquYSdgDgnP+Ec76Vcz7KOf9fchyzEIMOEyYXw6smsEuzKdfbUgDIEvYqC20mxFTHXMEdcprzer+X5pcRiqdwYNCGLnGGZTlR+2VxUv1WmYdOaNUq2EzaioRdI9oicm2gLolVp+uxmLZ1WXBJBmFv5gZgEioVg7NN3xQR+1PnF6BWsbyb4krAoFU13eZpXRlyCnnM2dHh6Wk/em3GjDivB7l6so+7Q7DoNXDlRHtDDsFCys26ef3ayvi3LtETLUvYRfGSW9gBwY7xVOCxS3cncqU8ClWn64uWt3ZZcH0xXPWHqJn7xGTTadE3hcf+1IUFvGHQDlueDXMlYNSpm27ztK5kj0iTOHHdh/0DtqqOm8+Kee6Su+CGZyEmPEGMdK7d0Bx2mhGOp9aknh27vgSHWYeBDtNKxF5GFsOl+SD0GlWmI6OcONt0mXmfpQhEkuizG9Gm12DGL4+wL4Xj6DCv7yK9tasNnANjC9VtoHqCcUX0K3FZDBs+3XHaF8H52YDi0hyzMWhJ2IsylDMibT4QxbQvgv0D1e2km3RqaFQsU6SUSKXxn755FH/xswsVHWfCHcJong1NKTMm1445ft2H/QN2MMbQldnsKh1hXZpfxubONqhrUMjhbNMXzLnPxR9JoN2oxaZ2g4wR+9p2AuWyVcyMqcZnl/rEKCFid1n0G17Yn76wAACKS3PMpt7j8ZpO2F0WPUw6dUYgpRzwff3VReyMsUyREgBcnFtGNJHGyxNepMtM4wvFkgU3NIczdxorwr4YiuOKJ4QDg8LarUYNDFpVmVZMsCY2DCBZMeXnsVuNWvTYjJj1y3PLX40NMthhgk6jylhV6yEQTYgtDZp78xQQrBhvKF6yV1Ejee3KInraDatqP5SGQatGpMnSHeuK1C9Eyow5PrkErZphV4+16mNnNwI7NeUHIGzMXigzfU662GRnxEj02AzQqhmuZGXGHMvy1wHhtXVZDSX7xfgjCcwFojUUdh2WY8mSt46pNMdyLIl2oxY9Nnki9mQqDX8kse6IXaNWYbOrLROxJ1Jp/OzMbEWeuxKKkySkXHa5UlFrwbQvgkGHWdZ6jI2GUatGlDZPizPkMGUi3+PXfdjZ0y7LlPp2ozaTFXNqyge9RvjzvDThLev3xzMZMWuFXaMW/PBsz/65S24YtepV+wNdFkPJiH1sQdo4Xfs8clBuLrtUddpu1KKn3QhPMF61jyhZYdXYIFu72nBpbhmeYAy/8aVX8JF/PYb//dPzZf++EtoJSLjaNn5bgVlfBD226gfFbGSMOrJiSjLkNGNyMYxYMoVTUz7sr9KGkcjuyX5yyo+bhjsw5DDhpfHyZpeOu0NgYvOvfAw7VoqrOOd49tICbtvsgF6zclHqtOpLCrvUmraWVgyAkj57ILoi7JvED+ZclXZMsQZg5bK124IZfxRv/8KLODHpw83DHfjmy9dwWrwLK4WShj10Wjd2kVIylcZcIIoem6HRS6kpBq2KNk9LMeQwIZHieObCAqKJdNUZMRLtorBH4ilcml/G3j4bDo068crE4qq8+UJc84bQ024sePcgVc2m0xwTnhAmFyO4c9vqTADJiinWjOzS/DKMWjV6axTlOMUeI6V8dn92xC5+MKvNjJEjWt4mXvAYgO9/9Fb8ywcOwmHW4zOPnC6r7cFixopRhscObNx+MQvLMaQ5FB+xG2jztDRSM7AfHp8GsOJRV4sk7Odm/UilOfb0tePQqAPLsSTOzgRK/v6sL1pUbIecZkQTacwvR/HsRaHD5V05BRndVgMiiRSWi7RavbywjC1dbTVrbSo1AiuV8igJu1W0YgBgpsoipSWps2MV+cx3bnXhzx7chUd/93bs7m2H1aDFZx7YgZNTfny7jGlWUi/23NF8zYhzg1sx0r7MpnZlR+xGSncsjdSm9pkLbjjbdOizy3O1bzfpEIgmMpk2e/ttuGWkAwBwZLy0zz4biGBTkVtKKTPmiieEZy8uYNRlXpOHnmncVMSOuVTDjBigfCsmO2LvFj+Ys1VuoBZrAFYuGrUKv3loKPM6AODBfT24ZaQDf/mzCyX3DjzBOCx6zSqLrFnRaVSwm7QbtkhpRrTuanX3uVEwatVUeVqKToseBq0K8VQa+/rtsu2mtxu14FwQ8S6rHl1WAzotBmzpbCu5gZpOc8z5o9jUXixiF0T8/OwyXrmyiLu2rS3IWGkrkF98fOE43Muxmm2cAsJto0WvKRnlZQu7QauGs01XtRUjeew2k7zRMmMMf/bgbgSiSXzv9amij232Ida5dFo27lDrTMSudGEXN0/lmPdQDk0p7IyxTAWqXP46sFJ9+tK4F3v6Vo57aNSB164sIp4s7LN7QjEkUrzoLWVPuxE6jQrfefU64sk07tq2ti9GqX4xUkvaLTWM2AFh47BUZJst7ACwqd1YtRWzGErArFPLkuWUy5YuCwYdJhwvMVdWKe0EJDqt+g0r7LO+CKwGTVV9npoBg1aNNAcSKRL2omSEXaaMGGBlzmIkkcLevpVBDreOOhBJCBk452cD+B8/PodHRH9fQsoGKSbsKhXDYIcJlxeCMGrVuGm4Y81jpM2uQiPyLtWwR0w2QvVpqXTHJHRqVWamoxy57Evh2kbL+/ttOHbdVzRy8obiitg4lXC16eHeoLNPp31RxW+cAvVv3du0wr65sw1aNcMeGYU9ezhydsR+87ADjAH/6Zuv4/6/fQFfOXwFDz8/sep3pUi11JtUai1w66gjr4dr1mtg0WsKpqddnl9Gm16DnhpvNpXTVsAvVp1KVpgQsUequt1cDMXzTk6Si/0DdriXYxlvNxfOOaaWwuhuV46w99iMmF+OlZXZVW9mWiCHHcgetkHCXpT/+MYRfO8jt8p6C5fdHXJPVsRuN+tw11YX7GYdPvPADrzrQC+u5XRqnBW95e4Sgitt/OazYSS62gsXKV0Ue8TUukrPaSltxQjtBFb+/r02I0LxFALR9Q9PXgrHazq1SLLuCtkxk4sRLEeT2NWzdvRes9LfYUQqzWVr+SAns/6I4jNiAMCoE6S2XhuoTSvs7UYt9soYrUvHBIQCo9z2oV/9rZvw5B/cid++YwR7etsRiqdWDTCY80eh06hKlqHv6rFCq2Z4U5FOdl1FipQuzwczedq1xNmmhy+cKNpjRGoAJiFlBM1WsYFa64h9e7cVeo0qk/mUy5kZoYhJjhYVG4U+u7BpP7m0dtCLRCCawGOnZuq2uQcIIrcUTrRWxC7DoPVyaFphrwWSSGXbMPmQ7JTsSU4z/ig2tRtKRtJv39ODFz95d+bDlg+hrcDaaHlsIQhvKI7dvbUXHSlVUCrWyUeusEsf0Gp89sVQbSN2nUaFG3rbC0bsZ2f8UKtYzfcw6omUDjy1VPi8/MkjZ/A7/3YcJybzX/BqgZRBpfSqUwDQSx47Rez1x6BV48F9PXjXgd6ij8v0hM/q+zLrK++WUqVimcyXQnRaDVhYjq6Jnh4/NQvGgLfs6i75PNVSTmHLGmEXUz2n15kZE02kEI6nap6Rsn/AhjMzAcTyRE9npgPY0tlWk6ycRrGp3QgVKyzsz1xcwI9OCPPnpRa69UAKAHqKpAgrBSNtnjaWv33PfrwpT355Nr3icOrsiH22RA57JXRZ9Uik+Jpo+fHTM7hxqKPkhUEOXJbSs08D0dXC7rLooVGxdUfsclSdlsP+ATviyTTOz67u2sk5x9kZv6L8dUC4S+m2GjC1uNaKCcWS+MwPz2DUZcb+ARuePF++sD99YR7BIhXSpZgtM+FACdDmaROgVavQZzdmGnql0hzzgahsm0D5ipQuzS/j0nwQb9uzSZbnKIWU7lcoMyad5gjkROxqFcMmmwHTRW75c/n0D07hY/92DED9uioW2kBdWI7BE6yP1VVv+uymvBH75564hGlfBP/nl/fgvl3dOD8bKOvCPB+I4kNfO4p/fGZs3Wua9kXAGOoSqDQao06yYuqTmUTCvk6ye8J7gjEk01y26jnpjZ692SXZMPftrr0NA2Q1AisQsQfjSaQ5YDWsrhDts5lWzaMtxYlJPx4/NYvXry1hKVR9y95y2NRuRLfVsGYD9Wxm41RZETsA9HUYMZWzeXpm2o+vHr6CX795ADcOdeCeHcKdajl2jPTef/z07Lo3XGf9Ebja9NBplC9DFLE3CcMOobc653ylLFqmyGN7twUuix7/9xcXERXLkB8/PYubhzvQaalPdCNUf6oKdniU+tbnDhDvtRsritilnu6ff/ISFsNSxF775lv7B2w4Prk6Yj87LTR627FJORunEn12E2YD0VXV0z8+OQO1iuGT920HIAyIGegwlSXsk6Ktc80bLqtBXj5mWqQ4CQD0YhEfeewbnEGHGcuxJBZD8ZWqU5l29816Df76V/fi0nwQf/mzi7g0H8TYQhAP7OmR5fjlwBhbVX3KOceFuZUPcHZnx2x6bUbML0eLtl/IJhBJwGLQ4IXLHvzi7ByA2nvsgCDsk4uRVZvDZ2b8GHaaYTE0f1fHXPrtRnC+OhX17EwA27otmYszYwx3b+/E4TFPyeyN64thMAZoVAyPnZpd15pm/JGWyIgBKGJvGqSGXle94UwVo5y7+3dudeEDhwbxlcNX8Nkfn4WKAffVIRsmm+zq068evor7Pv8CXhfH+QUihSN2zssbuCGN1vv1mwbgbNPhMdFuyj1mLZCGn2en952dCWCngvLXs8nksi8Kws45x5kZP3ZtWm073bOjE7FkGofHig+XmVwKo9tqwG2bnXj8dOX579KdbitkxABZLQUo3XFjI/WEv+YNYdYXgV6jkr0j4afu34HNnW04Mu7FLSMOuCz1LXOXIvbJxTD+6ucXAQBHry4CWD09KZs+8dZ6yle4GEZCujh0WQ34yJ2jmeNp1LV/W97Q2w6DVoXvvnYdnHP4wnFMLUWwW4H+OpCdyy6clxl/FL5wYs1G8c3DDph1ajxVwo6ZWoyg327CA3s2YXIxgtPT5U2nkvCFE4gm0orv6iihVaugVTOyYjY6fXYhN/iqN4zZgOAVyl3mb9Sp8fl374NBq8IvH+iT9djl4LLo4AnG8cc/PA0VEwZwnJwSItxMZ0fT2ogdQFk+e/bF4X03D8LZpq9bV0WDVo0/ePNWPHl+AT8+NYtzok+spIrTbDa1G6BWsUxmzBlRiHfmXMh0GhXeuNWFpy/MwxOMYTmav/p4cimMvg4jfmlnN7RqhscrtGOkDfbeFrFiAMCgUSOaqE9WjLJ7ZdYQvUaNHpsxE7HXqt/F7t52HPuTN2c8unoiRewvXI7hfzy4C69eWcQx0YrJeOyG1W+hTe1GMIayMmOy2/4adWr8/a/vx3IVfWYq5cO3j+Dx03P400fP4lcPChdOpQq7Rq3CpnZDJtPq7EwAKpZ/o/ieHV346Zk5HPyfTwIAdGoVHvnYbRmbKpZMYS4QRb/dhHaTFrdvduKxU7P41P3byw5uZjPdUFsjYgcAQx0HWlcVsTPGfpUxdpYxlmaMHZRrUc3CsNOMq54QZv3Rks2/qsGk09S86Vc+pL43bxi04zduHsS+fhtm/FEsBKLwRxJQq9iaJmw6jQqdFn1ZEXvuBuwtIw68eWeXzK+iMGoVw1/9yh4sRxN4+PkJbGo3wNGmnK6OufTZjZmI/dyMHyOuNph0a2O7d+ztweffvQ9/9uAu/NEvbUM8lcYrV1YGzcz4ouAcmelfD+zpwbQvgpNlDgsXjiG1E2gdYa/neLxqrZgzAN4F4HkZ1tJ0DDpMmPCEsLAcU+Qm0A197ehpN+AvfvkGqFQM+8SmaycmfULLXkP+C06vzVhxxN4otnZZ8Lt3bwHnyo3WJfrtpozHfmY6gN0FXq9Oo8I79/fiNw8N4T/fNQpnmy5jVQErqY79ou325p1dYAx4/pK77LXM+CPQqUs3zVMS9RyPV5UVwzk/D6Ah0eRGYMhhzlgHcqU6biTeMNiBI5++J/P17t52qFUMJ6d88EeSBQW5127CqanSzaQCEeFv10hhB4CP3jWK87MBvPWG+lT1Noo+uwnzgRhmfBHMBaJlFWIxxrBjkxXnZrOEXbw4SBF7u1ELm1Fb0cDsGV8Um2yGmg1k34gYdGrldXdkjD3EGDvKGDvqdpd/Zd/ISJkxgPKnrAPChuP2bgtOTvrXtBPIptdmxKwvinS6eArcRojYASFj4Yu/8Qa8fW/96gQaQX+HEGE/cW4eQPl3KDt7rLg8H8zUJkwuRqBVr25m12HWwRsqX9gn3EEMdBTucKpEDBrVxkl3ZIw9yRg7k+ffg5U8Eef8Yc75Qc75QZer8JCJZmLIsfLGbJVNoH39Npyc9MEXjq8pTpLotRsRT6VX9avPhz+SgFbNMqP1iNoi5bL/7IxQCFZu64RdPe2Ip9IYdwvzdicXw+izm6DOirYdbXp4S0zckogn07g8H1RszUAhjLoN5LFzzu/lnO/O8+9H9VjgRqa/wwTJhVKix56Pvf02LMeSOD+3XDDSzuSyl9hAldr+tqqVV2+kXPZXry6iz25ck6paiJ2bBAGWfPbJpXDmWBIOsw7eIr37sxl3BxFPpTPHbRWM2ibJiml1DFo1etqNMGrVq0bEKRlpeHg8mS4asQOlUx4D4sxUoj50WQ3QqhlSaV5RIdaw0wyDVpXx2ScXwxl/XcLRpis6lCUbpdcMFKJphJ0x9h8YY1MADgF4nDH2c3mW1TwMOU3osZWenKQURlxtmRTHYh47ULpIKbefO1Fb1CqWSS+sRFTVKobt3VacmwkgGEtiKZxAf84EsA6zHkvhOFIl9lUA4NxsAAatCsPOtspeQJNj0DVJgRLn/IcAfijTWpqST9+/A6Eqhg00G2oVw56+dhwZ9xYUZbNeA5tJi+kSbQX8kUTdKk0JgT67Ede8Yezurax1ws4eKx47OYPrXikjZq0Vw7kwLMVZohbg/GwA27osqzz6VsCgUSO6UTZPieLs7m3HzSOORi+jrkhDxItF2702Y1kee24/d6K2SJF2pTbIzk1WBKLJTKFSbsTuaBMu0KXsGM45zs0qt9laMYw6Vd2smNYwhglZ2VemsF/Jmgmbj9yZqUTtkXL1OyucHSAJ8c/F1sq5Hrt05+UJxooOAp8Vm4+12sYpIHjsyTRHIpWGtsaN7kjYiYq5Y4sTHzg0iENF7lR67Ua8OOYB5zzv/gPna0frEbXnjVtdeOPWytONt3dbwBjw6pVFtOk1sOdk1Ej2S6mIXdo4bcWI3ZA10LrWwk5WDFExJp0Gn31wN+xF/PE+uwnheAo+cdJSLsGYMFqPhL05MOk0GHaakeaCT597sZYi9lxhT6bSqzZUz4uZNdu6W0/Ypbmn9chlJ2EnakImM6ZAyuNKAzC6aWwWJPsk14YBhKlXjK0dfv7uh1/GH/2/k5mvz80GMOQwrWke1woYNKKw12GgNQk7URNWBjsUF3aK2JsHyT7J3TgFhGwpu0mHxay2ApxznJn24wfHpzP931t14xRYidjrsYFKwk7UhFIRu9QAjAqUmoeViD1/lXWHWbeqrYA7GENM7C/zf39xEcvRBK55wy25cQqszD0lYSeaFptJC6NWjdkSVgxF7M3DGwbteONWF+7Ykn/zNbetgHS3dvNwB5656Ma/vnwdALCjRYW9nnNPSdiJmsAYQ3e7AXOB/EOtA5kJTCTszYLFoMU3PnQTNnfmrxjNbSsgCfsfv3UHXBY9PveEMDe31a2YerTuJWEnakanRY+FQP4Oj4VmphLNi2DFrJxvaajH5s42/N7dm5FIcdhNWnRXmEOvFHb1WHHiv70Zd2x21vy5SNiJmlE0Yo8moGJAW57RbERz4jDr4YskkBSHX08tRdBh1sGs1+DdNw5goMOEvf22lumrlItWrYLNpIOmxjnsABUoETWky2rAfCCat0jJL3Z2bKUJOkrH0Sb1i0nAZdFjaimSyY7SaVT43kcPQaOiWLIe0F+ZqBldVgNiyXTGdsmG+sQoj9wipamcvu2dFgM1fasTJOxEzeiyCmXm83l8duoTozwcZuF8e0MxcM4xvRTJTG0i6gsJO1EzpE2yfD479YlRHlKHR28wnslhz520RNQHEnaiZkjDjufzCDtF7MrDkWXFSKmOJOyNgYSdqBmdkhXjzyfsSeoTozBsYr8YbzCWJexkxTQC+mQRNUOvUcNu0mJ+ebWwSy17qZ2AspD6xXhDcRjEHHaptQRRX0jYiZrSZTVgzr968zSaSCOeSpMVo0AcZqH6lAOwm7Qwt2AXx40A/dWJmtJlNWAhJ2IPRKlPjFKRGoGF4imyYRoIeexETem2GjCX47H7qU+MYnG26eENxdbksBP1hYSdqCldVj08wVimzBygzo5KpsOsgycYF3PYSdgbBQk7UVO62g1I89WTdfxhEnal0mHWwR9JiDnsZMU0ChJ2oqZ0WdbmspPHrlycbSstAyhibxwk7ERN6W5fW326Mu+UhF1pdIhtBQDKYW8kJOxETZGKlBbyCbuBkrKURnaTr16K2BsGCTtRU5xmPTQqtiZib9Nr6tKXmqgvkhVjN2nRRjnsDaOqTxZj7K8YYxcYY6cYYz9kjNnkWhihDFQqhk6LflWRUiCSJH9doUgRO9kwjaXakOkJALs553sAXALw6eqXRCiNzpwiJX8kAQvZMIrEZtJBxWjjtNFUJeyc819wzpPily8D6Kt+SYTSyC1Sopa9ykWtYrh52IFDo45GL6WlkTNs+hCA7xb6IWPsIQAPAcDAwICMT0tsdLqsehwZ9wAA0mmOGX8Eu1p0Un0r8O2Hbmn0ElqekhE7Y+xJxtiZPP8ezHrMfwWQBPCtQsfhnD/MOT/IOT/ocrnkWT3RFHS1GxCIJhGJp/DE+XlMLUXwlp3djV4WQSiWkhE75/zeYj9njH0AwNsA3MM553ItjFAOUpHSXCCKv3vqMoYcJjy4r6fBqyII5VJtVsx9AD4J4B2c87A8SyKUhlSk9JV82G4AAAUKSURBVG+vXMPZmQA+9qbNlOpIEDWk2k/X3wOwAHiCMXaCMfZPMqyJUBjSUOuvHr6K/g4j3rm/t8ErIghlU9XmKed8s1wLIZSLNPs0meb4nTdthpaidYKoKZRMTNScNr0GJp0aHWYd3nWAMmIJotaQsBM1hzGGT/zSNmzttlC0ThB1gISdqAsfvG240UsgiJaBwieCIAiFQcJOEAShMEjYCYIgFAYJO0EQhMIgYScIglAYJOwEQRAKg4SdIAhCYZCwEwRBKAzWiE67jDE3gGvr/HUnAI+My2kWWvF1t+JrBlrzdbfiawYqf92DnPOSAy0aIuzVwBg7yjk/2Oh11JtWfN2t+JqB1nzdrfiagdq9brJiCIIgFAYJO0EQhMJoRmF/uNELaBCt+Lpb8TUDrfm6W/E1AzV63U3nsRMEQRDFacaInSAIgihCUwk7Y+w+xthFxtgYY+xTjV5PLWCM9TPGnmGMnWeMnWWMfVz8fgdj7AnG2GXxv/ZGr1VuGGNqxthxxthj4tfDjLFXxNf8XcaYrtFrlBvGmI0x9j3G2AXxnB9S+rlmjP0X8b19hjH2bcaYQYnnmjH2FcbYAmPsTNb38p5bJvB3oradYowdqOa5m0bYGWNqAP8A4H4AOwG8lzG2s7GrqglJAH/IOd8B4BYAHxNf56cAPMU53wLgKfFrpfFxAOezvv4LAH8jvuYlAB9uyKpqy98C+BnnfDuAvRBev2LPNWOsF8DvATjIOd8NQA3gPVDmuf4agPtyvlfo3N4PYIv47yEAX6zmiZtG2AHcBGCMcz7BOY8D+A6ABxu8JtnhnM9yzo+J/78M4YPeC+G1fl182NcBvLMxK6wNjLE+AA8A+JL4NQNwN4DviQ9R4mu2AngjgC8DAOc8zjn3QeHnGsLkNiNjTAPABGAWCjzXnPPnASzmfLvQuX0QwDe4wMsAbIyxTet97mYS9l4Ak1lfT4nfUyyMsSEA+wG8AqCLcz4LCOIPoLNxK6sJnwfwCQBp8WsHAB/nPCl+rcTzPQLADeCrogX1JcaYGQo+15zzaQB/DeA6BEH3A3gdyj/XEoXOraz61kzCzvJ8T7EpPYyxNgDfB/D7nPNAo9dTSxhjbwOwwDl/PfvbeR6qtPOtAXAAwBc55/sBhKAg2yUfoqf8IIBhAD0AzBBsiFyUdq5LIev7vZmEfQpAf9bXfQBmGrSWmsIY00IQ9W9xzn8gfnteujUT/7vQqPXVgNsAvIMxdhWCxXY3hAjeJt6uA8o831MApjjnr4hffw+C0Cv5XN8L4Arn3M05TwD4AYBbofxzLVHo3Mqqb80k7K8B2CLunusgbLg82uA1yY7oLX8ZwHnO+eeyfvQogA+I//8BAD+q99pqBef805zzPs75EITz+jTn/H0AngHwK+LDFPWaAYBzPgdgkjG2TfzWPQDOQcHnGoIFcwtjzCS+16XXrOhznUWhc/sogPeL2TG3APBLls264Jw3zT8AbwVwCcA4gP/a6PXU6DXeDuEW7BSAE+K/t0LwnJ8CcFn8b0ej11qj138XgMfE/x8B8CqAMQD/D4C+0eurwevdB+CoeL4fAWBX+rkG8FkAFwCcAfBNAHolnmsA34awj5CAEJF/uNC5hWDF/IOobachZA2t+7mp8pQgCEJhNJMVQxAEQZQBCTtBEITCIGEnCIJQGCTsBEEQCoOEnSAIQmGQsBMEQSgMEnaCIAiFQcJOEAShMP4/gUzuthkpuqoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1089e8f28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "array([[ 10.],\n",
       "       [  1.]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim_length = 100\n",
    "x = np.array([[10., 1]]).T\n",
    "P = np.diag([5., 1.])\n",
    "F = np.array([[1., 1],\n",
    "             [0., 1]])\n",
    "R = np.array([[5., 0], [0, 5.]])\n",
    "Q = np.array([[.01, 0.0],\n",
    "              [0., .01]])\n",
    "H = np.array([[1., 0], [0., 1.]])\n",
    "\n",
    "xs, m = generate_simulation([10., 1], process_var=.01, sensor_var=5., length=sim_length)\n",
    "\n",
    "fx, fP = kalman(x, P, m, R, Q, F, H)\n",
    "\n",
    "plt.plot(range(sim_length), xs[:, 0], label=\"calculation\")\n",
    "plt.plot(range(sim_length), m[:, 0] ,'g', label=\"measurement\")\n",
    "plt.ylabel(\"distance\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.legend(loc='lower right', shadow=False)\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "\n",
    "plt.plot(range(sim_length), xs[:, 0], label=\"calculation\")\n",
    "plt.plot(range(sim_length), fx[:, 0], 'r--', label=\"filter\")\n",
    "plt.ylabel(\"distance\")\n",
    "plt.xlabel(\"time\")\n",
    "plt.legend(loc='lower right', shadow=False)\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "\n",
    "plt.plot(range(sim_length), fP[:, 0, 0])\n",
    "plt.show()\n",
    "\n",
    "# print(fx[:, 0])\n",
    "# print(xs - fx[:, 0])\n",
    "plt.plot(range(sim_length), xs[:, 0].flatten() - fx[:, 0].flatten())\n",
    "plt.show()\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kalman Filter Improvement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Major problem of the Kalman filter is that its performance depends a lot on proper estimation of process and noise covariance matrices Q & R. In practice these parameters are hard to estimate correctly, so filter has sub-optimal performance.\n",
    "\n",
    "In order to fix that, one of the proposals is to use neural network to correct kalman filter estimation.\n",
    "\n",
    "Network accepts $K$ and $y$ of the filter as the input and produces state correction $X_{mod}$, which is used to adjust state calculated by the filter.\n",
    "\n",
    "This idea is based on an article \"Kalman Filtering Compensated by Radial Basis Function Neural Network for Seam Tracking of Laser Welding\" (link below), which uses that technique to successfully improve Kalman Filter for Laser Welding process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kalman_training(x, P, measures, R, Q, F, H, B=0, u=0):\n",
    "    xs, cov = [], []\n",
    "    ys, Ks = [], []\n",
    "    for z in measures:\n",
    "        # predict\n",
    "        x = np.dot(F, x)\n",
    "        P = np.dot(F, P).dot(F.T) + Q\n",
    "\n",
    "        #update\n",
    "        S = np.dot(H, P).dot(H.T) + R\n",
    "        K = np.dot(P, H.T).dot(inv(S))\n",
    "        Ks.append(K.dot(np.array([[1], [0]])))\n",
    "        y = z - np.dot(H, x)\n",
    "        \n",
    "        ys.append(y)\n",
    "        x += np.dot(K, y)\n",
    "        P = P - np.dot(K, H).dot(P)\n",
    "        \n",
    "        xs.append(x)\n",
    "        cov.append(P)\n",
    "    return np.array(xs), np.array(cov), np.array(ys), np.array(Ks)\n",
    "\n",
    "\n",
    "sim_length = 1000\n",
    "x = np.array([[10., 1]]).T\n",
    "P = np.diag([5., 1.])\n",
    "F = np.array([[1., 1],\n",
    "             [0., 1]])\n",
    "R = np.array([[5., 0], [0, 5.]])\n",
    "Q = np.array([[.01, 0.0],\n",
    "              [0., .01]])\n",
    "H = np.array([[1., 0], [0., 1.]])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating Test Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run simulation l times for each sample until we reach a sample size specified below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 28%|██▊       | 142/500 [00:10<00:25, 14.20it/s]/Users/nick/anaconda3/lib/python3.6/site-packages/tqdm/_monitor.py:89: TqdmSynchronisationWarning: Set changed size during iteration (see https://github.com/tqdm/tqdm/issues/481)\n",
      "  TqdmSynchronisationWarning)\n",
      "100%|██████████| 500/500 [00:37<00:00, 13.21it/s]\n"
     ]
    }
   ],
   "source": [
    "import tqdm\n",
    "l = 1000  # Sample length\n",
    "N = 500 # Sample amount\n",
    "h1 = 50 # Number of neurons\n",
    "\n",
    "output_corr = []\n",
    "input_set = []\n",
    "\n",
    "for j in tqdm.tqdm(range(N)):\n",
    "    sim_length = 1000\n",
    "    x = np.array([[10., 1]]).T\n",
    "    P = np.diag([5., 1.])\n",
    "    F = np.array([[1., 1],\n",
    "                 [0., 1]])\n",
    "    R = np.array([[3., 0], [0, 3.]])\n",
    "    Q = np.array([[.01, 0.0],\n",
    "                  [0., .01]])\n",
    "    H = np.array([[1., 0], [0., 1.]])\n",
    "    xs, m = generate_simulation([10., 1], process_var=.01, sensor_var=3., length=sim_length)\n",
    "\n",
    "    fx, fP, Zs, Ks = kalman_training(x, P, m, R, Q, F, H)\n",
    "\n",
    "    correction = fx - xs\n",
    "    output_corr.append(correction.squeeze())\n",
    "    \n",
    "    inp = np.column_stack((Zs.squeeze(), Ks.squeeze()))\n",
    "    inp = np.array([inp for i in range(h1)])\n",
    "    input_set.append(inp)\n",
    "\n",
    "\n",
    "output_corr = np.array(output_corr)\n",
    "input_set = np.array(input_set)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural Network Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we use information from the article to create our neural network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = tf.placeholder(tf.float64, (N, h1, l, 4), name=\"x\")\n",
    "y = tf.placeholder(tf.float64, (N, l, 2), name=\"y\")\n",
    "W = tf.Variable(initial_value=np.random.randn(1, h1, l, 4))\n",
    "b = tf.Variable(initial_value=np.random.randn(1, l, 1))\n",
    "n = tf.exp(-tf.square(tf.norm(W - x, axis=1) * b))\n",
    "n = tf.reshape(n, (N, l, 2, 2))\n",
    "W2 = tf.Variable(initial_value=np.random.randn(1, l, 2))\n",
    "\n",
    "b2 = tf.Variable(initial_value=np.random.randn(1))\n",
    "out1 = tf.reduce_sum(W2*n[:,:,:,0], axis=2)\n",
    "out2 = tf.reduce_sum(W2*n[:,:,:,1], axis=2)\n",
    "\n",
    "out = tf.stack([out1, out2], axis = 2) + b2\n",
    "loss = tf.norm(y - out, ord=1)\n",
    "\n",
    "sgd = tf.train.RMSPropOptimizer(0.01)   # Learning rate\n",
    "step = sgd.minimize(loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train the Network"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We train the network trying to minimize loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: 591035.653106\n",
      "1: 549731.882457\n",
      "2: 522906.728771\n",
      "3: 503582.957777\n",
      "4: 488748.033239\n",
      "5: 476967.481459\n",
      "6: 467174.784179\n",
      "7: 458819.972324\n",
      "8: 451535.443752\n",
      "9: 445035.547035\n",
      "10: 439126.532735\n",
      "11: 433739.984486\n",
      "12: 428818.486771\n",
      "13: 424286.656795\n",
      "14: 420084.377217\n",
      "15: 416194.248795\n",
      "16: 412611.938753\n",
      "17: 409301.927858\n",
      "18: 406266.895894\n",
      "19: 403471.498664\n",
      "20: 400924.678517\n",
      "21: 398605.56388\n",
      "22: 396511.620445\n",
      "23: 394627.693255\n",
      "24: 392950.740841\n",
      "25: 391451.555363\n",
      "26: 390128.13236\n",
      "27: 388956.01068\n",
      "28: 387930.507074\n",
      "29: 387038.09779\n",
      "30: 386272.819012\n",
      "31: 385618.032477\n",
      "32: 385060.307688\n",
      "33: 384594.952936\n",
      "34: 384206.563997\n",
      "35: 383888.055959\n",
      "36: 383626.334888\n",
      "37: 383416.973692\n",
      "38: 383246.15842\n",
      "39: 383111.846082\n",
      "40: 383005.369447\n",
      "41: 382918.568282\n",
      "42: 382851.005882\n",
      "43: 382804.110335\n",
      "44: 382758.811399\n",
      "45: 382724.679157\n",
      "46: 382693.452271\n",
      "47: 382670.361672\n",
      "48: 382652.66863\n",
      "49: 382636.53439\n",
      "50: 382624.751275\n",
      "51: 382613.051332\n",
      "52: 382604.840244\n",
      "53: 382596.349349\n",
      "54: 382595.166754\n",
      "55: 382583.348514\n",
      "56: 382575.424466\n",
      "57: 382569.281113\n",
      "58: 382565.635516\n",
      "59: 382562.895458\n",
      "60: 382562.222466\n",
      "61: 382556.465232\n",
      "62: 382553.217866\n",
      "63: 382545.213314\n",
      "64: 382543.864991\n",
      "65: 382543.441842\n",
      "66: 382539.336273\n",
      "67: 382531.37871\n",
      "68: 382529.781138\n",
      "69: 382525.325108\n",
      "70: 382531.508406\n",
      "71: 382520.571109\n",
      "72: 382516.456812\n",
      "73: 382510.461991\n",
      "74: 382502.884607\n",
      "75: 382518.145669\n",
      "76: 382490.7672\n",
      "77: 382477.419248\n",
      "78: 382467.375863\n",
      "79: 382460.239525\n",
      "80: 382455.860687\n",
      "81: 382449.357781\n",
      "82: 382448.790913\n",
      "83: 382447.35864\n",
      "84: 382451.05326\n",
      "85: 382485.995776\n",
      "86: 382514.546691\n",
      "87: 382549.827621\n",
      "88: 382548.339554\n",
      "89: 382507.152689\n",
      "90: 382465.362833\n",
      "91: 382428.737159\n",
      "92: 382419.144128\n",
      "93: 382399.687197\n",
      "94: 382411.851506\n",
      "95: 382387.832376\n",
      "96: 382388.010365\n",
      "97: 382393.228162\n",
      "98: 382397.48685\n",
      "99: 382404.384453\n",
      "100: 382411.199157\n",
      "101: 382404.020454\n",
      "102: 382390.837367\n",
      "103: 382382.731821\n",
      "104: 382360.068794\n",
      "105: 382347.859035\n",
      "106: 382341.26983\n",
      "107: 382336.568017\n",
      "108: 382333.075572\n",
      "109: 382327.935202\n",
      "110: 382330.816332\n",
      "111: 382333.11875\n",
      "112: 382329.976503\n",
      "113: 382325.865091\n",
      "114: 382322.568585\n",
      "115: 382308.450471\n",
      "116: 382314.230855\n",
      "117: 382290.855233\n",
      "118: 382285.09952\n",
      "119: 382279.876721\n",
      "120: 382278.224157\n",
      "121: 382280.89171\n",
      "122: 382271.349085\n",
      "123: 382263.536145\n",
      "124: 382260.062678\n",
      "125: 382264.176503\n",
      "126: 382256.038034\n",
      "127: 382248.26538\n",
      "128: 382246.893751\n",
      "129: 382243.992964\n",
      "130: 382234.754598\n",
      "131: 382223.147479\n",
      "132: 382218.551219\n",
      "133: 382213.253738\n",
      "134: 382213.038272\n",
      "135: 382214.515745\n",
      "136: 382206.400395\n",
      "137: 382236.740597\n",
      "138: 382207.81773\n",
      "139: 382187.424038\n",
      "140: 382178.652118\n",
      "141: 382169.676907\n",
      "142: 382166.20523\n",
      "143: 382157.355972\n",
      "144: 382160.736543\n",
      "145: 382154.811498\n",
      "146: 382149.460455\n",
      "147: 382150.057361\n",
      "148: 382151.464492\n",
      "149: 382145.754233\n",
      "150: 382135.686529\n",
      "151: 382129.591159\n",
      "152: 382119.251569\n",
      "153: 382114.595736\n",
      "154: 382119.592577\n",
      "155: 382102.434775\n",
      "156: 382093.66704\n",
      "157: 382096.885888\n",
      "158: 382105.610146\n",
      "159: 382097.900069\n",
      "160: 382103.380748\n",
      "161: 382091.300417\n",
      "162: 382091.800671\n",
      "163: 382084.444299\n",
      "164: 382068.769279\n",
      "165: 382053.710913\n",
      "166: 382046.902333\n",
      "167: 382042.638862\n",
      "168: 382034.322343\n",
      "169: 382038.695617\n",
      "170: 382035.425092\n",
      "171: 382052.8115\n",
      "172: 382049.125267\n",
      "173: 382051.608191\n",
      "174: 382030.397615\n",
      "175: 382022.783557\n",
      "176: 382018.914196\n",
      "177: 382002.544222\n",
      "178: 382005.954998\n",
      "179: 381995.958645\n",
      "180: 381997.553895\n",
      "181: 381991.933887\n",
      "182: 381984.560113\n",
      "183: 381990.201118\n",
      "184: 381984.462746\n",
      "185: 381979.964019\n",
      "186: 381971.027858\n",
      "187: 381960.805161\n",
      "188: 381958.716584\n",
      "189: 381953.814897\n",
      "190: 381951.71403\n",
      "191: 381951.711882\n",
      "192: 381939.102549\n",
      "193: 381945.17858\n",
      "194: 381933.624457\n",
      "195: 381931.653722\n",
      "196: 381938.282396\n",
      "197: 381940.019818\n",
      "198: 381931.789425\n",
      "199: 381925.424545\n",
      "200: 381914.427227\n",
      "201: 381906.029093\n",
      "202: 381902.943346\n",
      "203: 381913.708518\n",
      "204: 381898.073208\n",
      "205: 381899.405445\n",
      "206: 381898.992211\n",
      "207: 381889.948706\n",
      "208: 381882.029224\n",
      "209: 381875.459418\n",
      "210: 381876.374791\n",
      "211: 381866.692125\n",
      "212: 381862.153518\n",
      "213: 381861.111647\n",
      "214: 381860.661104\n",
      "215: 381866.772094\n",
      "216: 381860.955877\n",
      "217: 381852.883215\n",
      "218: 381857.516842\n",
      "219: 381847.867101\n",
      "220: 381839.824308\n",
      "221: 381823.604425\n",
      "222: 381825.982531\n",
      "223: 381835.661002\n",
      "224: 381826.883725\n",
      "225: 381817.268849\n",
      "226: 381818.179834\n",
      "227: 381823.444444\n",
      "228: 381815.420371\n",
      "229: 381832.963937\n",
      "230: 381823.903723\n",
      "231: 381812.83508\n",
      "232: 381809.684461\n",
      "233: 381805.022844\n",
      "234: 381794.753439\n",
      "235: 381795.442983\n",
      "236: 381790.336691\n",
      "237: 381783.778985\n",
      "238: 381779.283941\n",
      "239: 381780.550854\n",
      "240: 381787.783051\n",
      "241: 381782.362658\n",
      "242: 381786.220705\n",
      "243: 381784.479778\n",
      "244: 381779.749335\n",
      "245: 381767.549582\n",
      "246: 381752.370548\n",
      "247: 381764.754639\n",
      "248: 381745.618196\n",
      "249: 381745.014372\n",
      "250: 381738.029888\n",
      "251: 381734.924746\n",
      "252: 381747.934119\n",
      "253: 381752.718813\n",
      "254: 381762.226027\n",
      "255: 381758.49846\n",
      "256: 381741.283697\n",
      "257: 381743.032349\n",
      "258: 381746.52556\n",
      "259: 381727.714062\n",
      "260: 381713.981466\n",
      "261: 381714.882546\n"
     ]
    }
   ],
   "source": [
    "epochs = 1000\n",
    "sess = tf.InteractiveSession()\n",
    "init = tf.global_variables_initializer()\n",
    "sess.run([init])\n",
    "\n",
    "for i in range(epochs):\n",
    "    loss_, _ = sess.run([loss, step], feed_dict={x: input_set, y: output_corr})\n",
    "    print(str(i) + \": \" + str(loss_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you run the training code above, you will notice that the loss is quite huge and not at all acceptable. This means that the network has failed to train on a given test data. \n",
    "\n",
    "I was not able to repeat experiment result of the article, which is due to several issues.\n",
    "First of all, neural network was poorly described in the article, which makes it hard to infer its configuration. I lack experience in this domain (machine learning), so to analyze what went wrong and where is hard as it requires a lot of expert knowledge.\n",
    "\n",
    "Secondly, paper describes concrete process and concrete solution to predict that process more accurately. I do not see why the algorithm proposed cannot be used in general sense, but again, lack of expert knowledge leaves ma unsure.\n",
    "\n",
    "Quality of the paper leaves more to be desired. Notation is inconsitent and there are a lot of grammatical mistakes which makes it hard to understand what the authors meant. Also authors leave a lot of details out of the paper, such as the size of the training dataset, initial parameters, etc.\n",
    "\n",
    "Given more time and knowledge, this research has a great potential to enhance classic filtering algorithms, however I was not able to repeat the success of the authors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://ieeexplore.ieee.org/document/6353188/ - paper the research is based on"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
